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/*MDL_BEGIN 

MDL_SOURCE_F I LE_NAME : R I CKS_R 9 6 MATH . H 

PKG_ACRONYM : R96MA! MDL_SOURCE_VERSION# : VI . 0 

MDLJDESCRIPTION: 
This module is an include file that provides some of the 
math support for TTI's REAL 9 6 math datatype. 
********************************************* 

* Copyright (c) 2001 by * 

* Touch Technologies Inc., San Diego, California. * 

* ALL RIGHTS RESERVED. * 
**************************************************************************** 

R96 datatype example implementation 
Data stored in Little endian 
» bits 0-63 are the Integer Part (IP) (signed, so 63 bits of precision) 

1 bits 64-95 are the Fractional Part (FP) (signed, modulus 100,000,000 - ie. 8 digits) 

! Sign is generally stored in both IP and FP, however, 

! determine the sign of a r96, as follows: 

1 If IP != 0 then IP is 2s complement and can be used to determine the sign 

1 If IP ===== 0 then FP is 2s complement 

! Because zero has no sign, to determine the sign if IP=0, you have to 

! look at the sign stored in fractional part. 

i 

/* In specific embodiments, routines are provided for unsigned mutliprecision 

y= 1 Integer Math (64, 96, and 128 bit) 

L~ MACROS and routines are organized herein so that 64, 96, and 12 8 bit variants are in order for each 
y routine or MACRO. This makes maintenance easier because the variants are very similar. 

Q Most math is done in unsigned 96bit integer format (scaled) . Conversion to 96bit unsigned integer state is 

done by removing the sign and then multiplying the IP by 100,000,000 and adding in the FP. For multiply 
™ and divide, some 128 bit math may be required to prevent bit loss. Optimizations have been performed to 
%y minimize the number of bits required, especially when multiplication is required. 

yg There are MACROS to perform 196 compares, additions and subtractions. Most MACROS assume that the 
parameters are of type (unsigned long *) 



Miscellaneous MACROS 



I There are MACROS to perform comparisons such as ==, !=, >, <, and compares against ZERO. There are also 
MACROS to copy an integer and to negate an integer. 



Bit shift MACROS 



= 3 : 



There are MACROS to perform right and left bit shift with single bit shifts and counted shifts. Counted 
y shifts MUST be 31 bits or less. There are MACROS to check for overflow as well used in multiplication. 



H* Addition MACROS: 



fy If a 2 operand ADD is needed then the SUM and the second ADDEND must be the same variable and NOT the 

first ADDEND, i.e. a = a + b is coded as add (b, a, a) . 
Zr\ There are MACROS to handle overflow as well necessary for multiply. The i##add32 variant adds a 32 bit 
I y integer to the appropriate extended precision integer providing a shortcut for that case 



Subtraction MACROs: 



These MACROs perform unsigned subtraction. 



Multiplication Routines: 



! The basic mechanism is to determine whether the multiplicand or the 

! multiplier is smaller and use that number as the multiplier. The 

I multiplier is then shifted left until bit ZERO is set. Each time 

! the multiplier is shifted left, the multiplicand is shifted right. 

1 When bit ZERO is set in the multiplier, the current shifted 

! multiplicand is added to the product. When the multiplier becomes 

! ZERO, the multiplication stops and the product is now complete. 

I We do the multiplication by BIG "switch" statement that does 8 

! bits at a time to save time over doing a single bit at a time. 

! The "ixxmultiply" and n ixxmultiply_oiverf low" routines are created by 

! a program named "create_ixxmultiply_r96math. c" . It does the 64, 96, 

! and 128 bit routines. 

1 

1 There are also routines to test for and handle overflow. Overflow is 

• identical to non-overflow except the shift and add routines called 

I check for the overflow condition. 
! 

! There are also some special MACROS that multiply specifically by 10, 

! 10,000, and 100,000,000. These speed up the various R96 conversions 

! and operations. 
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I NOTE: These routines are 3 operand routines ! I i One CANNOT use the same variable for the I 
product as either the multiplicand or the multiplier!!! 

! There is a set of routines to locate the highest set bit. These are 

! used to optimize overflow handling in multiplication. it turns out that 

! overflow can be easily predicted by adding the highest set bits and 

! comparing them against the number of the sign bit. >= will be overflow. 

! The final result is also compared to see if the sign bit is set which 

! implies overflow before sign handling is done. 



J Integer Divides 

! The algorithm shifts the divisor over until it is >= the dividend and 

i then shifts it back one bit. The divisor is then subtracted from the 

! dividend. The quotient gets the value of the number of bit shifts that 

! were done added to it. The difference of the subtract becomes the 

! remainder. This process is repeated until the remainder is less than 

! the divisor. 



R96 miscellaneous routines and MACROS 

! r96abs is a MACRO that converts an R96 number into a positive value 

! r96negate is a MACRO that negates an R96 number 

! r96copy is a MACRO that copies an R96 number from one location to another 

! r96cvtr96toi96 is a routine that converts an R96 number into an unsigned 

I 96 bit integer. The sign is returned as the function value (0 if positive 

I and 1 if negative) . The IP is multiplied by 100,000,000 and then the 

„ , I FP is added in. 

f* ! r96cvtr96toi96__4digitfp is a special case routine of r96cvtr96toi96 that 

O 1 converts an R96 number into an unsigned 96 bit integer when modulus 10,000 

^ ! of the FP == ZERO. The sign is returned as the function value (0 if 

! positive and 1 if negative). The IP is multiplied by 10,000 and then the 

yy ! FP is added in. 

J3 1 r96cvti96tor96 is a routine that converts an unsigned 96 bit integer 

^ ! into an R96 number. The sign is a parameter. The 196 value is 

yj ! divided by 100,000,000. The quotient becomes the IP and the remainder is 

fl I the FP. 

HI ! r96cvtasctobin converts an ASCIZ string into an R96 number. 

! r96cvtbintoasc converts an R96 number into an ASCIZ string. 
s ! I f this routine is to be used as part of the TTI R96 math stuff, 

Q 1 it- should have an optimization that uses a table lookup where the 

n 1 the binary is divided by 10,000 converting 4 digits at a time, 

i instead of 10 which only converts 1 digit at a time. 

H= l ! R96 division ~~ " 

f|j ! Divide is accomplished by converting the dividend and divisor to 96bit 

p ! integers and then doing a 96bit divide. The remainder is then converted 

^ ! to a 128 bit integer and multiplied by 100,000,000. This value is then 

fy ! divided by the original divisor using the least required bits to get the 

! FP portion of the result quotient. Overflow is tested by checking bit 

1 63 of the result quotient. Overflow can happen when dividing by a 

! fractional entity, i.e. 100 / .1 == 100 * 10. 

1 

! We call r96divide_nohilong as an optimization when the IP part of both 

! numbers is < 2^32. This routine is able to save 32 bits of precision 

! in the operations, 
i 

! Overflow handling is done by checking the unsigned final quotient 

• (IP portion) against 0x7FFFFFFFFFFFFFFF before signs are applied. 

! If the value is larger, then runtime$error is called to signal the 

! error. The runtime$error in the R96 math .H files is intended 

! to be over-written. The definition is conditional. 



R96 multiplication ~~ 

! The multiplication starts with a conversion to 128 bit unsigned 

! integers, a 128 bit multiply is then done. The product is 

! then divided by {100,000,000 * 100,000,000) to get the IP part 

• of the product. The remainder of that divide is divided by 
1 100,000,000 to get the FP part of the product. 

! We start by checking for possible optimizations. If both FPs == ZERO 

! then we call r96multiply_nofp which performs 164 math. If either 

• the multiplicand or the multiplier == ZERO, then we call 
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! r96multiply_nofp_multiplicand or r96multiply_nofp_multiplier as 

! appropriate. Finally, if modulus 10,000 of both FPs « 0 then we call 

! r96multiply_4digitfp. The last 3 optimization routines perform 196 bit 

! math instead of 1128 bit thus running faster. 

! 

I Overflow handling is done by checking each left shift and addition for 
! overflow as well as checking the unsigned final product (IP portion) 
! against 0x7FFFFFFFFFFFFFFF before signs are applied. If the value is 
! larger, then runtime$error is called to signal the error and processing 
I continues until a solution is arrived at. The runtime$error in the R96 
1 math .H files is intended to be over-written. The definition is 
! conditional. THIS IMPLIES THAT runtime$error CAN BE CALLED MULTIPLE 
I TIMES. Noticing an overflow does not directly stop the processing of 
I the multiplication. 
*/ 

# include tt ixxmath_nomultiply .h" 
#include vv ixxmath_multiply. h" 
# include u r96math.h" 
/*MDL_END*/ 



* COPYRIGHT * 

* Copyright (C) 2001 * 

* by Touch Technologies Inc., San Diego, California * 

* This software is furnished under a license and may be used and * 

* copied only in accordance with the terms of such license and with * 

* the inclusion of the above copyright notice. This software or any * 

* other copies thereof may not be provided or otherwise made * 
□ * available to any other person. No title or ownership of the * 

* software is hereby transferred * 
******************************************************^ 

0 FACILITY: SheerPower Language 
f! ABSTRACT: 

2 Tni s module provides the real integer 18.8 math class for 

y both Alpha & Intel Pentium Processors.*/ 

1 Module for the real integer 18.8 math class for Alpha & Intel Pentium 

#include "c-universal .h" 
^ #include <stdio.h> 
.J #include <string.h> 
sj # include M c~intmsg.h" 
r 7" # include <time.h> 
** #include <math.h> 
y # include "ricks_r96math.h" 
=» typedef struct 

Ti ^ f nt wholedi 9 its ? /* number of whole digits eg. 99.678 has 2 whole digits*/ 

y int fractdigits; /* number of fractional digits eg. 99.678000 has 3 fractional digits*/ 

int neg; /* if the number is negative*/ 

int rlen; /* result length eg. 99.678 has an rlen of 6*/ 

} realinfo; 
long ran_random () ; 
int randomize (); 
void scaledown (void *_a, int n) ; 
int ator(real *a, char *b) ; 

void rtoa (real *, void *, int fdigits, realinfo *rinfo, int flags); 
inline void scaleup (void *_a) / 
static real scale = 1; 

#define LARGE S T_NUMB E R ( ~ { ( (int64 ) 1 ) «63 ) ) 

/ / ' 1 ' 1 f ' t f 1 ' 1 ' ' f ' 1 1 L 1 ' ' ' ' L L > > ' > 1 1 f t t f f t 1 / / , t i t t t t t t , { 1 t f f t , t , i , , i , / 



MTHSDEXP DOUBLE EXPONENTIAL 

' ; ' 1 . 1 ' ' 1 f ' 1 ' ' ( ' 1 ' ! 1 1 ! 7 1 < / ' ' * 1 1 1 / / / n 1 1 / / / / / / / / / / / / / / / 7 / / / / / / / / / / / 

// Brief description: Return the exponential of the input value. 

// Expected: d = the address of a 'double' 

// Result: return = the exponential of 'd' as a double 

'a / U ' ' 'J- ' ' ' J^LJ 1 i ' 1 1 ' ' 1 ' 1 ' ' ' ' /'//////////////////////////////// / 

double routine mth$dexp (double *d) 

{ return exp (*d); // call exp() passing the value of 'd'} 

/ 1 ' 1 ' 1 1 ' 1 ! 1 ' 1 ' ' L 1 L 1 i ' L LU ' > ' 1 ' > > 1 ' 1 < J r i t t i f i t t i f < t , t , t , 1 , t f / , t } t i j , , 

n MTH$DLOG2 DOUBLE BASE 2 LOGARITHM 

' t / 1 ' > ' J L 1 ( ' '. f ' ' ' ' / / ' ' ' ' ' ' 1 ' 1 ' ' ' ' ' ' 1 1 f ' 1 / / ' > < >> * > ~ 1 / / / / / / / / / / / / / / 7 / 
// Brief description: Return the Base 2 logarithm of the input value. 

// Expected: d = the address of a 'double' 

// Result: return = the base 2 logarithm of 'd» as a double 
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t i t f i t t it tit / // // // // / / // // // // // // // // // // // // // // // // // // // // // // 
double routine mth$dlog2 (double *a) 
{ return log (*a) * 1.442695099;} 

/ 1 f L f i L L i I i i t f i i i f i / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / 



MTH$DNINT DOUBLE TO INTEGRAL PART 



/////////////////////////////////////////////////////////////////// 

// Brief description: Remove the fractional part from a double and return result. 

// Method: Convert the double to a large 64-bit integral value. This removes 

// the fractional part. Then convert the large 64-bit integer back 

// to a double, and return that value. 

// Expected: a = the address of a 'double' 

// Result: return = the integral part of 'a' as a double. 

'////////////////////////////////////////////////////////////////// 
double routine mth$dnint (double *a) 

{ int64 i = (int64) * a; // convert to a large interger 

double d = (double) i; // then to a double 

return d; // return integral part as a double} 

/ * 1 1 1 ' L f 1 1 f 1 t i i i * t i * i j i f / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / 



OTS$POWDD RAISE A DOUBLE TO A DOUBLE EXPONENT 

I 1 1 / * 1 1 / / i * t i ' * i t i t / / / / / / / / / / / / / / / / / / 7 i / i t / / / / / / / / / / / / / / / / 7 / / / / 7 i t i 

II Brief description: Raise a double 'x' to a double exponent 'y' . 
// Expected: x = the base double 

// : y = the exponent double 

// Result: return = double 'x' raised by 'y' 

ffttittitiittit i / / / i if tii ii / // // // // / // // // // // // // // // // // // // // // 

double routine ots$powdd (double x, double y) 
{ double r = pow(x, y) 
return r; } 

1 ' f f ' 1 1 1 ' i t * i ' ' i f j t i f t i > I I I t i I I I I I I I I I I i t I l f I t I I t i I i I ! I i ! j f j j t I t i ! i 



% n OTS$POWJJ RAISE AN INTEGER TO AN INTEGER EXPONENT 



^ 1 i f i t t i 1 t i i t / t t i i i i i i i i i i i i i i i it iii i i i i i f i i i i i i i i i t i i i i i t i i t ( i i i i j i 
■W // Bri ef description: Raise an integer 'x' to a integer exponent 'y' . 
» // Expected: 



% II 
*M II Result 



x = the base integer 

y = the exponent integer 

return = double 'x' raised by 'y' 

tttiitititftftitti t t i / / / / / // // // // // // // // // // // // // // // // // // // // / 
double routine ots$powjj (int x, int _y) 

P { 

double v = v; 

U| double x = _x; // convert to double's first for pow() 

double r = pow(x, y) ; 
s return r; } 

jF«& i 1 1 / 1 ' 1 ' 1 ' ' 1 1 1 1 L i ' t i i i ' i i f i i f f ' i t i i i i i t t i i i ( i i i i i i i i j / / / / / / / / / / / / / / / 
OTSSP OWDJ RAISE A DOUBLE TO AN INTEGER EXPONENT 

—7 — ; — : — r~-s — : — ; — : — : — r— ; — : — : — . — ■. — ■. — — -. — ; — ■ — ■ — — — . 



fx! ttiti/itftttiii/ii///// / t i t tit i i i i i i i i i t i i i f t t / i i i i i i i i i i i i i i { t i i f i 

If Brief description: Raise a double 'x' to an integer exponent 'y». 



// Expected 
// 

i // Result 



x = the base double 

y = the exponent integer 

return = double 'x' raised by 'y* 



b—s // ^c&uit; icLUiu = aouDie 'x- raisea oy 'y 

■WSF tllfllfl/itlilltilffllfi i ////// I / // // // // // // // // // // // // // // // // // 

ffj double routine ots$powdj (double x, int _y) 

{ double y = _y; // convert to double first 

double r = pow(x, y) ; 
return r; } 

y * ' ' 1 f ' ' f 1 1 ' ' f * f 1 f t I t f l * t I I I I I I j j j t { f I t t l I t l I t I I f l I j j j j j f f I l j l I l j / } i 

1 u MTHSDSQRT SQUARE ROOT OF A DOUBLE 

' ' ' ' ' ' 1 ' ' ' ' 1 / 1 ' 11 f 'I * I > < I l t I I I I i I I t i I l { / t t i { i i i t f , i f , i i , , ! i i , , t , ! , , 

II Brief description: The square root of a double. 

// Expected: d = the address of the double 

// Result: return = the square root of 'd' as a double 

1 1 / / ' ' 1 i i 1 1 * i * i / / f i t ' 1 i i / / t f / ! t i I I I I t I t I I / I i i i i i i i i i i t i i t i i t j i t i / / / 

double routine mth$dsqrt (double *d) 
{ return sqrt (*d) ; } 

t I ! I t t t I I I t I i t I / / / / / / / / / f ? I I I ! I I I I I I i f I t I I I t f I i t i j I t I I f I l l f l l j f i i f / 



MTH$DLOG NATURAL LOGARITHM 

1 ' '.' ' ' ' 1 1 ' ' 1 1 1 t 1 * f > ' ' ' t 7 7 / 7 t t t i i t t i i i t i i i i i i i / / / / / / / / / / 7 7 t t 7 7 , ~, 7 7 , , 7 

// Brief description: The natural logarithm of a double as a double. 

// Expected: d = the address of the double 

// Result: return = the natural logaithm of 'd' as a double 

1 1 1 ( 1 ' 1 1 1 1 1 ' 1 i I I I I I I I I I i l I I I I l I i I i 1 I I i i i i i I I i I I i I i i I i l I i ! I t i ( , i l l ! , 

double routine mth$dlog (double *d) 
{ return log (*d) ; } 

1 / 1 1 / 1 ' { 1 1 1 1 1 1 ' 1 L 1 ' t f 1 1 f t t 1 t 1 1 t 1 1 1 1 i t t t f 1 1 1 t / / 1 , / / / / / / / / / f _ , i f t f f _ f t , 
u MTH$DLOG10 COMMON LOGARITHM 



' { ' 1 1 1 1 1 1 1 ' ' 11 > ' > I > > 1 ' ' 1 ' l < * l 1 l f I 1 t i I t f i I I f t f i I t i l t l t 1 f 1 t 1 i i 1 ; , t , , / 

// Brief description: The common logarithm of a double. 

// Expected: d = the address of the double 

// Result: return = the common logarithm of 'd' as a double 

f ' ' ' f ' < ' ' ' f f ' * > ' > I i ' I I ' I l i t t f t I i I I t ( { t { i i i t f 1 t ( t t , t i t t t I I I I t I I I l l [ l 

double routine mth$dlogl0 (double *d) 
{ return loglO (*d> ; } 
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/////////////////////////////////////////////////////////////////// 



n MTH$DSIN THE SINE OF ITS RADIAN ARGUMENT 

/////////////////////////////////////////////////////////////////// 

// Brief description: return the sine of its radian argument 

// Expected: d = the address of the double 

// Result: return = the return the sine of its radian 

// argument ' d' as a double 

/////////////////////////////////////////////////////////////////// 
double routine mth$dsin (double *d) 

{ return sin (*d) ; // return the sine of its radian argument} 

/////////////////////////////////////////////////////////////////// 



MTH$DSINH THE HYPERBOLIC SINE OF ITS RADIAN ARGUMENT 



/////////////////////////////////////////////////////////////////// 

// Brief description: return the hyperbolic sine of its radian argument 

// Expected: d = the address of the double 

// Result: return = return the hyperbolic sine of its radian 

// argument 'd' as a double 

/////////////////////////////////////////////////////////////////// 
double routine mth$dsinh (double *d) 

{ return sinh (*d) ; // return the hyperbolic sine of its radian argument} 

/////////////////////////////////////////////////////////////////// 



MTH$DCOS THE COSINE OF ITS RADIAN ARGUMENT 



/////////////////////////////////////////////////////////////////// 
// Brief description: return the cosine of its radian argument 
// Expected: d = the address of the double 

// Result: return = return the cosine of its radian 

// argument ' d' as a double 

/////////////////////////////////////////////////////////////////// 
Ijl double routine mth$dcos (double *d) 

i"Z { return cos (*d) ; // return the cosine of the radian argument} 

ipS / // // // // // // // // // // // // // // // // // // // // // // // // // // // // // // // // // 

£ f n MTHSPCOSH THE HYPERBOLIC COSINE OF ITS RADIAN ARGUMENT 

: =ssr / // // // // // // // // // // // // // // // // // // // // // // // // // // // // // // // // // 

■„f\ II Brief description: return the hyperbolic cosine of its radian argument 
^ // Expected: d = the address of the double 

%g // Result: return = return the hyperbolic cosine of its radian 

yi // argument ' d' as a double 

:p /////////////////////////////////////////////////////////////////// 
O double routine mth$dcosh (double *d) 

ijft { return cosh (*d) ; // return the hyperbolic cosine of its radian argument} 

y s ////////////////////////////////////////////////////////////////// / 

s 1 u MTH$DTAN THE TANGENT OF ITS RADIAN ARGUMENT 

(StfB ; / // // // // // // // // //"> "/"/""/ //////////////////////////////////////////// 

1*3 // Brief description: return the tangent of its radian argument 

!;| // Expected: d = the address of the double 

IfT // Result: return = return the tagent of its radian 

^ // argument ' d' as a double 

S=j ■/////////////////////////////////////////////////////////////////// 

double routine mth$dtan (double *d) 
O { return tan (*d) ; // return the tangent of its radian argument} 

—7" / // // // // // // // // // // // // // // // // // // // // // // // // // // // // // // // // // 
li j n MTHSPTANH THE HYPERBOLIC TANGENT OF ITS RADIAN ARGUMENT 

/////////////////////////////////////////////////////////////////// 
II Brief description: return the hyperbolic tangent of its radian argument 
// Expected: d = the address of the double 

// Result: return = return the hyperbolic tangent of its radian 
// argument "d" as a double 

/////////////////////////////////////////////////////////////////// 
double routine mth$dtanh (double *d) 

{ return tanh (*d) ; // return the hyperbolic tangent of its radian argument} 

/////////////////////////////////////////////////////////////////// 



// MTH$DATAN THE ARC TANGENT OF ITS RADIAN ARGUMENT 



/////////////////////////////////////////////////////////////////// 

// Brief description: return the arc tangent of its radian argument 

// Expected: d = the address of the double 

// Result: return = return the arc tagent of its radian 

// argument 'd' as a double 

/////////////////////////////////////////////////////////////////// 
double routine mth$datan (double *d) 

{ return at an (*d) ; // return the arc tangent of its radian argument} 

/////////////////////////////////////////////////////////////////// 



MTH$DASIN THE ARC SINE OF ITS RADIAN ARGUMENT 



/////////////////////////////////////////////////////////////////// 

II Brief description: return the arc sine of its radian argument 

// Expected: d = the address of the double 

// Result: return = return the arc sine of its radian 

// argument 'd' as a double 

/////////////////////////////////////////////////////////////////// 
double routine mth$dasin (double *d) 

{ return asin (*d) ; // return the arc sine of its radian argument } 

/////////////////////////////////////////////////////////////////// 



Appendix, Page 5 of 43 



QIPLG Attorney Docket No:500.0028-30us 



n MTH$DACOS THE ARC COSINE OF ITS RADIAN ARGUMENT 



/////////////////////////////////////////////////////////////////// 

// Brief description: return the arc cosine of its radian argument 

// Expected: d - the address of the double 

// Result: return = return the arc cosine of its radian 

// argument 'd' as a double 

/////////////////////////////////////////////////////////////////// 
double routine mth$dacos (double *d) 

{ return acos (*d) ; // return the arc cosine of its radian argument} 

/////////////////////////////////////////////////////////////////// 



MTH$DABS THE ABSOLUTE VALUE OF A COMPLEX NUMBER 



/////////////////////////////////////////////////////////////////// 

// Brief description: return the absolute value of a complex number 

// Expected: d = the address of the double 

// Result: return = return the absolute value of a complex 

// number as a double 

/////////////////////////////////////////////////////////////////// 
double routine mth$dabs (double *d) 
{ return fabs (*d) ; } 

/////////////////////////////////////////////////////////////////// 



// LIB$EDIV EXTENDED PRECISION DIVIDE 

/////////////////////////////////////////////////////////////////// 
II Brief description: divide a 64 -bit quadword by an integer giving an integer quotent 
// and remainder. 

// Expected: a = integer divisor 

// b - address of 64-bit quadword 

// c address of integer quotent 

// d = address of integer remainder 

// Result: status from lib$ediv 

/////////////////////////////////////////////////////////////////// 
c_routine int routine lib$ediv (void *a, void *b, void *c, void *d) 
{ *(int *) c = (int) (*(int64 *) b / Mint *) a); 

Mint *) d = (int) (*(int64 *) b % *{int *) a); 

return 1 ; } 

ft L i L-LJ. i f i i t i L i i t f i t i / / / / / / iiii / / / / / // // // / / / /////// / / / // / / // / / / / / / 
// EMUL EXTENDED PRECISION MULTIPLY 

/'>'>//////////////////////////////////////////////,/////,/////,,// 
// Brief description: multiply two integers giving a quadword result. 
// Expected: 

// a integer multiplier 

// b integer multiplicand 

// c = integer addend 

// q = address of quadword result 

// Result: status from lib$emul 

/''//////////////////////////////////////////////////////////////// 
int routine lib$emul (int *_a, int *_b, int *c, void *_d) 
{ int64 a = *_a; 

int64 b = *Jb; 

int64 d = a * b; 

if (c) 

d += *c; 

*(int64 *) _d = d; 

return 1 ; } 

t t t t f t t f f t t t f t t f t t / t t ////// t ////////////////////////////////////// / 
n EDiV EXTENDED PRECISION DIVIDE 

///////////// i t f / t //)////// t ///// / } //////////////// / / // // // // // // / / 
// Brief description: divide a 64-bit quadword by an integer giving an integer quotent 
// and remainder. 

// Expected: a = integer divisor 

// b = address of 64 -bit quadword 

// c = address of integer quotent 

// d = address of integer remainder 

// Result: status from lib$ediv 

t 1 1 > t i * i > / t i t / i / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / i f / t / / / / / / / / / / / / / f t 
c_routine int routine ediv (int a, int *b, int *c, int *d) 
{ return lib$ediv (&a, b, c, d) ; } 

jiitttiiiiii i f i i / i i it f ii i i i i i i i i i i j t i i / i i i i i / i i j / i / i t { t ( f i t / / t i j { i i 
u EMUL EXTENDED PRECISION MULTIP 

/ / ' / 7 ' i ' * if ft t ' J''''/'/////////////// ///////////////////////// / i 

II Brief description: multiply two integers giving a quadword result. 
// Expected: 

// a integer multiplier 

// b = integer multiplicand 

// c = integer addend 

// q = address of quadword result 

// Result: status from lib$emul 

ttttfiitttit/tt t / / / / / / /t// / // // // / / // // // // // // // // // // // // // // // // 
c_routine int routine emul (int a, int b, int c, int *q) 



Appendix, Page 6 of 43 



QIPLG Attorney Docket No:500.0028-30us 



{ return lib$emul (&a, &b, &c, q) ; } 

/////////////////////////////////////////////////////////////////// 
// SUBX EXTENDED PRECISION SUBTRACT 

///// ////// //////////////////////////////////////////////////////// 
// Brief description: subtract two 64 -bit quadword giving quadword result 
// Expected: 

// a address of 64 -bit quadword 

// b s address of 64 -bit quadword 

// c address of 64 -bit quadword 

// Result: status from lib$subx 

/////////////////////✓///////////////////////////////////////////// 
c_routine int routine lib$subx (void *a, void *b, void *c) 
{ *(int64 *) c = *(int64 *) a - * (int64 *) b; 
return 1 ; } 

/////////////////////////////////////////////////////////////////// 
n LIB$ADDX EXTENDED PRECISION ADD 

/////////// //////////////////////////////////////////////////////// 
// Brief description: add two 64-bit quadword giving quadword result 
// Expected: 

// a = address of 64 -bit quadword 

// b address of 64 -bit quadword 

// c address of 64 -bit quadword 

// Result: status from lib$addx 

/////////////////////////////////////////////////////////////////// 
c_routine int routine lib$addx (void *a, void *b, void *c, void*) 
{ *(int64 *) c = *(int64 *) a + * (int64 *) b; 
return 1 ; } 



RANDOM NUMBERS 



I I I I I I I I I I I I I I i 1 I I t l I I 1 l I l I I I l I I I I I I I t i I I I I I I I I i ! I I I I I I t i I I I I I I I { t 1 

; MTH$RANDOM GET A RANDOM VALUE 

/////////////////////////////////////////////////////////////////// 

// Brief description: Get a random value modulus 'maxvalue' or if no 'maxvalue' 

// return a value between 0 and 1 . 

// Expected: maxvalue = address of an integer 

// Result : floating random number 

/////////////////////////////////////////////////////////////////// 

real routine mth$random (int *maxvalue) 

{ int r = ran_random {); /* get random integer*/ 

real one = 1; 

real f; 

if (*maxvalue) 

{ r %= *maxvalue; /* remainder after divide by maxvalue*/ 

f = (real) (r + 1) ; /* relative to one*/ } 

else 

{ f = (real) r; 

while (f >= one) /* make the random number less than one*/ 

f /= 10; } 
return f; /* return floa/real*/} 

#include <stdio.h> 
/* random.c: 



An improved random number generation package. 

/* In addition to the standard rand { ) /srand ( ) like interface, this package also has a special state info 
interface. The initstateO routine is called with a seed, an array of bytes, and a count of how many 
bytes are being passed in; this array is then initialized to contain information for random number 
generation with that much state information. Good sizes for the amount of state information are 32, 64, 
128, and 256 bytes. The state can be switched by calling the setstateO routine with the same array as 
was initiallized with initstateO. By default, the package runs with 128 bytes of state information and 
generates far better random numbers than a linear congruential generator. If the amount of state 
information is less than 32 bytes, a simple linear congruential R.N.G. is used. 

Internally, the state information is treated as an array of longs; the zeroeth element of the array is the 
type of R.N.G. being used (small integer) ; the remainder of the array is the state information for the 
R.N.G. Thus, 32 bytes of state information will give 7 longs worth of state information, which will allow 
a degree seven polynomial. (Note: the zeroeth word of state information also has some other information 
stored in it see setstateO for details). 

The random number generation technique is a linear feedback shift register approach, employing trinomials 
(since there are fewer terms to sum up that way) . In this approach, the least significant bit of all the 
numbers in the state table will act as a linear feedback shift register, and will have period 2^deg - 1 
(where deg is the degree of the polynomial being used, assuming that the polynomial is irreducible and 
primitive) . The higher order bits will have longer periods, since their values are also influenced by 
pseudo-random carries out of the lower bits. The total period of the generator is approximately 
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deg*(2**deg - 1); thus doubling the amount of state information has a vast influence on the period of the 
generator. 

Note: the deg*(2**deg - 1) is an approximation only good for large deg, when the period of the shift 
register is the dominant factor. With deg equal to seven, the period is actually much longer than the 
7* (2**7 - 1) predicted by this formula. 

For each of the currently supported random number generators, we have a break value on the amount of state 
information (you need at least this many bytes of state info to support this random number generator) , a 
degree for the polynomial (actually a trinomial) that the R.N.G. is based on, and the separation between 
the two lower order coefficients of the trinomial.*/ 

//////////////////////////////// 



/////// 


/////// 


/ / / / / 


/ / 


/////////////// 


#def ine 


TYPE_0 


0 


/* 


linear congruential*/ 


#def ine 


BREAK_0 


8 






#def ine 


DEG_0 


0 






#def ine 


SEP_0 


0 






ttdefine 


TYPE_1 


1 


/* 


x**7 + x**3 + 1*/ 


#def ine 


BREAK_1 


32 






ttdef ine 


DEG_1 


7 






ttdefine 


SEP_1 


3 






#def ine 


TYPE_2 


2 


/* 


x**15 + x + 1*/ 


#def ine 


BREAK 2 


64 






#def ine 


DEG 2 


15 






#define 


SEP_2 


1 






#def ine 


TYPE_3 


3 


/* 


x**31 + x**3 + 1*/ 


#def ine 


BREAK_3 


128 






ttdef ine 


DEG_3 


31 






ttdef ine 


SEP_3 


3 






f=k ttdef ine 


TYPE_4 


4 


/* 


x**63 + x + 1*/ 


r™* #def ine 


BREAK_4 


256 






ZZ ttdef ine 


DEG 4 


63 






U ttdefine 


SEP 4 


1 







/* * Array versions of the above information to make code run faster -- relies 

* on fact that TYPE_i == i.*/ 

ttdefine MAX_TYPES 5 /* max number of types above*/ 

static int degrees [MAX_TYPES 3 = { DEG_0 , DEGJL, DEG_2,DEG_3, DEG_4 } ; 
static int seps [MAX_TYPES] = { SEP_0, SEP_1, SEP_2,SEP_3, SEP_4}; 
/* * Initially, everything is set up as if from : 

* initstate( 1, fcrandtbl, 128 ); 

* Note that this initialization takes advantage of the fact that ran_srandom ( ) 

* advances the front and rear pointers 10*rand__deg times, and hence the 

* rear pointer which starts at 0 will also end up at zero; thus the zeroeth 

* element of the state information, which contains info about the current 

* position of the rear pointer is just 

* MAXJTYPES* (rptr - state) + TYPE_3 == TYPE_3.*/ 
static long randtbl [DEG_3 + 1] = { TYPE_3, 

0x9a319039, 0x32d9c024, 0x9b663182, 0x5dalf342, 

0x48f 340fb, 
0x946554fd, 
0x2d436b86, 
0x904f35f7, 
0xac94efdc, 
0x8a88d77b, 



0xde3b81eO, 0xdf0a6fb5, 
0x7449e56b, OxbebldbbO, 



0x8c2e680f , 0xeb3d799f , 
0xda672e2a, 0xl588ca88, 



0xd7158fd6, 
0x36413f93, 



0x6fa6f051, 
0xc622c298 , 



0xfl03bc02, 
0xab5c5918, 
0xbllee0b7, 
0xe369735d, 
0x616e6b96, 
Oxf 5a42ab8, 
0xf5ad9d0e, 0x8999220b, 0x27fb47b9}; 
/* * fptr and rptr are two pointers into the state info, a front and a rear 

* pointer. These two pointers are always rand_sep places aparts, as they cycle 

* cyclically through the state information. (Yes, this does mean we could get 

* away with just one pointer, but the code for randomO is more efficient this 

* way) . The pointers are left positioned as they would be from the call 

* initstate( 1, randtbl, 128 ) 

* (The position of the rear pointer, rptr, is really 0 (as explained above 

* in the initialization of randtbl) because the state table pointer is set 

* to point to randtbl [1] (as explained below).*/ 
static long *fptr = &randtbl [SEP_3 + 1] ; 

static long *rptr = &randtbl [1] ; 

/* * The following things are the pointer to the state information table, 

* the type of the current generator, the degree of the current polynomial 

* being used, and the separation between the two pointers. 

* Note that for efficiency of random () , we remember the first location of 

* the state information, not the zeroeth. Hence it is valid to access 

* state [-1], which is used to store the type of the R.N.G. 

* Also, we remember the last location, since this is more efficient than 
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* indexing every time to find the address of the last element to see if 

* the front and rear pointers have wrapped.*/ 
static long *state = &randtbl [1] ; 

static int rand_type = TYPE_3 ; 

static int rand_deg = DEG_3 ; 

static int rand_sep = SEP_3; 

static long *end_ptr = fcrandtbl [DEG_3 + 1] ; 

/* * ran_s random: 

* Initialize the random number generator based on the given seed. If the 

* type is the trivial no-state- information type, just remember the seed. 

* Otherwise, initializes stated based on the given "seed" via a linear 

* congruential generator. Then, the pointers are set to known locations 

* that are exactly rand_sep places apart. Lastly, it cycles the state 

* information a given number of times to get rid of any initial dependencies 

* introduced by the L.C.R.N.G. 

* Note that the initialization of randtbl [] for default usage relies on 

* values produced by this routine.*/ 
void ran_s random (unsigned x) 

{ register int i, j; 

if (rand_type TYPE_0) 

{ state [0] = x; } 

else 

{ j - l; 

state [03 = x; 

for (i = 1; i < rand_deg; i++) 
{ state [i] = 1103515245 * state [i - 1] + 12345/ } 

fptr = &state [rand_sep] ; 
a& rptr = &state[0]; 

for {i=0; i<10* rand_deg; 
!f ran_random () ? 

□ }} 

fjj /* * initstate: 

^ * Initialize the state information in the given array of n bytes for 

y * future random number generation. Based on the number of bytes we 

,0 * are given, and the break values for the different R.N. G. ' s, we choose 

«l * the best (largest) one we can and set things up for it. ran_srandom{ ) is 

* then called to initialize the state information. 

J I * Note that on return from ran_srandom() , we set state [-1] to be the type 
. * multiplexed with the current value of the rear pointer; this is so 
^ * successive calls to initstate () won't lose this information and will 
s=? * be able to restart with setstateO. 

lj * Note: the first thing we do is save the current state, if any, just like 
^ * setstateO so that it doesn't matter when initstate is called. 
7~ * Returns a pointer to the old state.*/ 

y char *ran_initstate (unsigned seed, /* seed for R. N. G. */ 
3 cnar *arg_state, /* pointer to state array*/ 

int n) /* # bytes of state info*/ 

^ { register char *ostate = (char *) (&state [-1] ) ; 
if (rand_type == TYPE_0) 

state C - 13 = rand_type; 
else 

state [-1] = MAX_TYPES * (rptr - state) + rand_type; 
if (n < BREAK_1 ) 

{ if (n < BREAK_0 ) 

{ fprintf (stderr, 

"initstate: not enough state (%d bytes) with which to do jack; ignored, 
n) ; 

return 0; } 
rand_type = TYPE_0; 
rand_deg = DEG_0; 
rand_sep = SEP_0 ; } 

else 

{ if (n < BREAK_2) 

{ rand_type = TYPE_1; 

rand_deg = DEG__1; 
rand_sep = SEP_1; } 

e] se 

{ if (n < BREAK_3 ) 

{ rand_type = TYPE_2; 

rand_deg = DEG_2 ; 
rand_sep = SEP_2; 
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} 

else 

{ if (n < BREAK_4 ) 

{ rand_type « TYPE_3; 

rand_deg = DEG_3 ; 
rand_sep - SEP_3; 

} 

else 

{ rand_type = TYPE_4; 

rand_deg = DEG_4; 

rand_sep = SEP_4; 
}}}} 

state = &(((long *) arg_state) [1] ) ; /* first location*/ 

end_ptr = &state [rand_deg] ; /* must set end_ptr before ran_srandoin*/ 

ran_srandom (seed) ; 

if (rand_type == TYPE_0) 

state [-1] = rand_type; 
else 

state [-1] = MAX__TYPES * (rptr - state) + rand_type; 
return (ostate) ; } 
/* * setstate: 

* Restore the state from the given state array. 

* Note: it is important that we also remember the locations of the pointers 

* in the current state information, and restore the locations of the pointers 

* from the old state information. This is done by multiplexing the pointer 

* location into the zeroeth word of the state information. 

* Note that due to the order in which things are done, it is OK to call 
|s& * setstate () with the same state as the current state. 

f*. * Returns a pointer to the old state information.*/ 

;Jj!]J char *ran_setstate (char *arg_state) 

W { register long *new_state = (long *) arg_state; 

J3 register int type = new_state[0] % MAX_TYPES; 

7=i register int rear = new_state[0] / MAXJTYPES; 

~a? char *ostate = (char *) (&state [-1] ) ; 

yl if (rand^type == TYPE_0) 

f s t ate [ - 1 ] = rand_type ; 

:*Z else 

U\ state [-1] = MAX_TYPES * (rptr - state) + rand_type; 

s switch (type) 

{ case TYPE_0: 

case TYPE_1 : 
yd case TYPE_2: 

\dk case TYPE_3: 

L s case TYPE_4 : 

rand_type = type; 
rand_deg = degrees [type] ; 
S j; rand_sep = seps [type] ; 

^ break; 

default: 

fprintf (stderr, 

"setstate: state info has been munged; not changed."); } 
state = &new_state [1] ; 
if (rand_type != TYPE_0) 

{ rptr = &state [rear] ; 

fptr = &state [ (rear + rand__sep) % rand_deg] ; } 
end_j?tr = Estate [rand_deg] ; /* set end_ptr too*/ 
return (ostate) ; } 
/* * random: 

* If we are using the trivial TYPE_0 R.N.G., just do the old linear 

* congruential bit. Otherwise, we do our fancy trinomial stuff, which is the 

* same in all ther other cases due to all the global variables that have been 

* set up. The basic operation is to add the number at the rear pointer into 

* the one at the front pointer. Then both pointers are advanced to the next 

* location cyclically in the table. The value returned is the sum generated, 

* reduced to 31 bits by throwing away the "least random" low bit. 

* Note: the code takes advantage of the fact that both the front and 

* rear pointers can't wrap on the same call by not testing the rear 

* pointer if the front one has wrapped. 

* Returns a 31-bit random number.*/ 
long routine ran_random ( ) 

{ long i; 
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if (rand_type == TYPE_0) 

{ i = state [0] = (state [0] * 1103515245 + 12345) & 0x7fffffff ; } 

else 

{ *fptr += *rptr; 

i = (*fptr >> 1) & 0x7fffffff ; /* chucking least random bit*/ 
if (++fptr >= end_j>tr) 
{ fptr = state; 

++rptr; } 

else 

{ if (++rptr >= end_j)tr) 

rptr = state; 

} } 

return ( i ) ; } 
static unsigned char const ascii_to_ebcdic [] = 
{ 0, 01, 02, 03, 067, 055, 056, 057, 

026, 05, 045, 013, 014, 015, 016, 017, 

020, 021, 022, 023, 074, 075, 062, 046, 

030, 031, 077, 047, 034, 035, 036, 037, 

0100, 0117, 0177, 0173, 0133, 0154, 0120, 0175, 

0115, 0135, 0134, 0116, 0153, 0140, 0113, 0141, 

0360, 0361, 0362, 0363, 0364, 0365, 0366, 0367, 

0370, 0371, 0172, 0136, 0114, 0176, 0156, 0157, 

0174, 0301, 0302, 0303, 0304, 0305, 0306, 0307, 

0310, 0311, 0321, 0322, 0323, 0324, 0325, 0326, 

0327, 0330, 0331, 0342, 0343, 0344, 0345, 0346, 

0347, 0350, 0351, 0112, 0340, 0132, 0137, 0155, 

0171, 0201, 0202, 0203, 0204, 0205, 0206, 0207, 
M 0210, 0211, 0221, 0222, 0223, 0224, 0225, 0226, 
^ 0227, 0230, 0231, 0242, 0243, 0244, 0245, 0246, 
^ 0247, 0250, 0251, 0300, 0152, 0320, 0241, 07, 
U 040, 041, 042, 043, 044, 025, 06, 027, 
,f| 050, 051, 052, 053, 054, Oil, 012, 033, 
,15 060, 061, 032, 063, 064, 065, 066, 010, 

070, 071, 072, 073, 04, 024, 076, 0341, 
yjl 0101, 0102, 0103, 0104, 0105, 0106, 0107, 0110, 
:S 0111, 0121, 0122, 0123, 0124, 0125, 0126, 0127, 
;~-:f 0130, 0131, 0142, 0143, 0144, 0145, 0146, 0147, 
yH 0150, 0151, 0160, 0161, 0162, 0163, 0164, 0165, 
B 0166, 0167, 0170, 0200, 0212, 0213, 0214, 0215, 

f ^ 0216, 0217, 0220, 0232, 0233, 0234, 0235, 0236, 
M 0237, 0240, 0252, 0253, 0254, 0255, 0256, 0257, 
Ly 0260, 0261, 0262, 0263, 0264, 0265, 0266, 0267, 
y. 0270, 0271, 0272, 0273, 0274, 0275, 0276, 0277, 
C] 0312, 0313, 0314, 0315, 0316, 0317, 0332, 0333, 
IU 0334, 0335, 0336, 0337, 0352, 0353, 0354, 0355, 
O 0356, 0357, 0372, 0373, 0374, 0375, 0376, 0377}; 
s=r ; static unsigned char const ebcdic_to_ascii [] = 
= { 0, 01, 02, 03, 0234, 011, 0206, 0177, 

0227, 0215, 0216, 013, 014, 015, 016, 017, 

020, 021, 022, 023, 0235, 0205, 010, 0207, 

030, 031, 0222, 0217, 034, 035, 036, 037, 

0200, 0201, 0202, 0203, 0204, 012, 027, 033, 

0210, 0211, 0212, 0213, 0214, 05, 06, 07, 

0220, 0221, 026, 0223, 0224, 0225, 0226, 04, 

0230, 0231, 0232, 0233, 024, 025, 0236, 032, 

040, 0240, 0241, 0242, 0243, 0244, 0245, 0246, 

0247, 0250, 0133, 056, 074, 050, 053, 041, 

046, 0251, 0252, 0253, 0254, 0255, 0256, 0257, 

0260, 0261, 0135, 044, 052, 051, 073, 0136, 

055, 057, 0262, 0263, 0264, 0265, 0266, 0267, 

0270, 0271, 0174, 054, 045, 0137, 076, 077, 

0272, 0273, 0274, 0275, 0276, 0277, 0300, 0301, 

0302, 0140, 072, 043, 0100, 047, 075, 042, 

0303, 0141, 0142, 0143, 0144, 0145, 0146, 0147, 
0150, 0151, 0304, 0305, 0306, 0307, 0310, 0311, 
0312, 0152, 0153, 0154, 0155, 0156, 0157, 0160, 
0161, 0162, 0313, 0314, 0315, 0316, 0317, 0320, 
0321, 0176, 0163, 0164, 0165, 0166, 0167, 0170, 
0171, 0172, 0322, 0323, 0324, 0325, 0326, 0327, 
0330, 0331, 0332, 0333, 0334, 0335, 0336, 0337, 
0340, 0341, 0342, 0343, 0344, 0345, 0346, 0347, 
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0173, 0101, 0102, 0103, 0104, 0105, 0106, 0107, 
0110, 0111, 0350, 0351, 0352, 0353, 0354, 0355, 
0175, 0112, 0113, 0114, 0115, 0116, 0117, 0120, 
0121, 0122, 0356, 0357, 0360, 0361, 0362, 0363, 
0134, 0237, 0123, 0124, 0125, 0126, 0127, 0130, 
0131, 0132, 0364, 0365, 0366, 0367, 0370, 0371, 
060, 061, 062, 063, 064, 065, 066, 067, 
070, 071, 0372, 0373, 0374, 0375, 0376, 0377}; 
/* * lib$tra_asc_ebc -- translate ASCII to EBCDIC 
** J 

int routine lib$tra_asc_ebc (DESC_S *source, DESC_S *dest) 
{ unsigned char *asc, *ebc; 

int si en, dlen; 

int i ; 

split$desc (source, &slen, &asc) ; 
split$desc (dest, &dlen, &ebc) ; 
for (i = 0; i < slen && i < dlen; i++) 
ebc[i] = ascii_to_ebcdic [asc [i] ] ; 

return SS$_N0RMAL;} 
/* * lib$tra_ebc_asc translate EBCDIC to ASCII 

■k-k J 

int routine lib$tra_ebc_asc (DESC_S *source, DESC_S *dest) 
{ unsigned char *asc, *ebc; 

int slen, dlen; 

int i; 

split$desc (source, &slen, &ebc) ; 

split$desc (dest, &dlen, &asc) ; 

for (i = 0; i < slen && i < dlen; i++) 

asc [i] = ebcdic_to__ascii [ebc [i] ] ; 
return SS$_N0RMAL; } 
/* * Convert packed decimal to leading separate 
* 

* string: pointer to output ascii string buffer 

* packed: pointer to input packed decimal string*/ 
int routine cvtps(char *string, char *packed) 

{ int leading = 1; 
int negative = 0; 
char buf [32+1] ; 
char *s = buf; 
int i, C; 

for (i = 0; i < 32; i++) { if (i&l) /* odd is lower*/ 

c = packed [i/2] & OxOf? 
else /* even*/ 

c = (packed [i/2] >>4) & OxOf; 
if (leading && c == 0) 

continue; /* skip leading 0's*/ 
leading =0; /* nolonger leading*/ 
if (c == 10 || c == 12 || c == 14 || c == 15) 

break; /* plus*/ 

else if (c == 11 | | c == 13) { negative = 1; 

break; } 
else if (c > 9) 

break; 
*s++ = (char) (c + 1 0 ' ) ; 
} 

*S++ = *\0' ; 

if (negative) 

*string++ = '-'; /* no unary '+'*/ 
for (i = 0; buf [i] ; i++) 

*string++ = buf [i] ; 
if (i == 0) 

*string++ = '0'; /* at least one zero*/ 
*String++ = 1 \0 1 ; 
return SUCCESS ; } 
/* * Convert leading separate to packed decimal 
**/ 

int routine cvtsp(char *packed, char *string, int len) 
{ int negative = 0; 

char buf [32+1] ; 

char *s = buf; 
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int i; 

for (i = 0; i < len s < buf +32; i++) { if {string [i] <= 1 1 ) 

continue; 
else if (string [i] == '.*) 

continue ; 
else if (string [i] == '-') 

negative = 1; 
else if (string [i] == ' + ') 

negative = 0; 
else if (string [i] >= '0' && string [i] <= '9') 

*s++ = string [i] ; 

} 

*s++ = ' \0' ; 

len = strlen(buf); 

/* * if first zero padding required*/ 

if ((len&l) == 0) { for (i = 0; buf [i] ; i++) ; 

for {; i >= 0; --i) 
buf [i + 1] = buf [i] ; 

buf[0] = ' 0'; 

len = st r len (buf) ; 

} 

for (i = 0; i < len; i++) { if (i&l) /* odd is lower*/ 

packed [i/2] |= (buf [i] &0x0f ) ; 
else 

packed[i/2] = (char) ( (buf [i] &0x0f) «4) ; 

} 

M> packed [i/2] |= negative ? 13 : 12; 

f-s return (i/2)+l; /* return length of packed decimal string*/} 

SI /* * Convert packed decimal to real 

&p * packed: pointer to input packed decimal string 

.A * result: pointer to output result real*/ 

H J int routine cvtpdr(real ♦result, char *packed, int digits) 

yO { char string [32 + 1] ; /* maximum packed decimal length*/ 

Q int i, n; 

rz /* * Convert packed decimal to leading separate 

m */ 

s cvtps (string, packed); 

/* * set packed decimal precision to 'digits' fractional digits*/ 

for (i = 0; string[i]; i++) ; 
|£J for (n = 0; n < digits && i >= 0; --i, n++) 

La string [i+1] = string [i] ; 

Ls string [i+l] = 

L!: /* * Convert ascii string to real*/ 

rj return ator (result, string);} 

/* * Convert real to packed decimal 

■ *"•••'" * 

* packed: pointer to output packed decimal string 

* digits: sizeof output buffer (as packed decimal digits) 

* val : pointer to real value to convert*/ 

int routine cvtrpd (char *packed, int digits, int precision, real *val) 
{ realinfo info; 

char buf [32+1]; /* maximum size of a packed decimal field*/ 

char pbuf [32+1] ; /* packed buffer*/ 

int p_digits; 

DESC_S desc; 

int i; 

set$desc (&desc, sizeof (buf), buf); 

/* * Convert real to ascii string*/ 

rtoa(val, &desc, precision, &info, 0) ; 

/* * Convert ascii string to packed decimal*/ 

p_digits = cvtsp(pbuf, buf, 0) ; 

digits - (digits/2) +1; /* convert digits to bytes + sign*/ 

/* * fill packed decimal string with required 

* leading packed zeros.*/ 
for (i = p_digits; i < digits; i++) 

*packed++ = ' \0 ' ; 
/* * now fill with converted packed number*/ 
for (i = 0; i < p_digits; i++) 

*packed++ = pbuf [i] ; 
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return SUCCESS;} 
/* * cvtlp -- convert long to packed decimal 
**/ 

int routine cvtlp (int *expr, int digits, char *buf) 
{ real r = *expr; 

return cvtrpd(buf, digits, 0, &r) ; } 
/* * cvtpl -- convert packed decimal to long 
* */ 

int routine cvtpl (int digits, char *buf, int *result) 
{ int stat; 
real r ; 

stat = cvtpdr(&r, buf, digits); 
* result « (int)r.ip; 
return stat; } 

DEBUG DUMP HEX BYTES 

/ / i / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / i i / / / / / / / i f f > i i i i i i f i 

a DHEXB DUMP HEX BYTES OF A MEMORY LOCATION 

/////////////////////////////////////////////////////'////////'//'' 

// Brief description: This function is a debug function. It will dump in hexadecimal 

// format with ascii characters at the right, a block a memory 

// starting at address 'a' for bytes 'n'. Unprintable character 

// hex values are displayed but their ascii value is ' . ' . 

// This function is called when debugging the math routines 

// at a very low level. It is called to dump the contents of 

// 'real' 64-bit integers {and their 128-bit intermediate values) 

// to check that the hex values are indeed correct. 

// Any block of memory may be dumped with this function; not 

// just 'real' 64-bit (or intermediate 128-bit) values. I have 

// used this function to dump 'pcode' blocks of memory when I 

// was not sure what was being generated. 

// Expected: a = the start address (any type) 

// n number of bytes to dump. 

// Result: return = void 

//////////////////////////////////////////////////////////''//''''' 

void routine dhexb (void *a, int n) 

{ unsigned char *s = (unsigned char *) a; 

char msg[256] ; 

char buf [32] ; 

int sn = n; 

int i ; 

for (; n > 0; n -= i) 

{ sprintf (msg, 11 %04d ", sn - n) ; 

for (i = 0; i < 16 && i < n; i++) 
{ sprintf (buf, "%02X ", s [i] ) ; 

s treat (msg, buf) ; } 
for (; i < 16; i++) 

s treat (msg, " " ) ? 
s treat (msg, " " ) ; 

for (i =0; i < 16 && i < n; i++, s++) 
{ sprintf (buf, »%c" , *s < ' • || *s > ? '.' : *s) ; 

strcat (msg, buf) ; } 
logf ( "%s\n" , msg) ; 

////^ )////////////////////////////////////////////////////////////// 
// Brief description: add 2 reals and place result in a third 
// Expected: 1st arg = unused 

// a the address of a quadword (integer pointer) 

// b the address of a quadword (integer pointer) 

// c the address of a quadword (integer pointer) 

// 

// Result: c = a + b; 

/////////////////////////////////////////////////////////////////// 
void routine addm (int, int *a, int *b, int *c) 
{ *{int64 *) c = *{int64 *) a + *(int64 *) b;} 



MORE MATH AND COMPARE FUNCTIONS 



/////////////////////////////////////////////////////////////////// 



a SUBM SUBTRACT 2 QUADWORDS AND PLACE RESULT IN A THIRD 

~ 7 7 1 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 T~~~7~~l 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 ? 7 7 7 7 7 7 7 7 7 1 7 7 7 7 7~ 

II Brief description: subtract 2 reals and place result in a third 
// Expected: 1st arg = unused 

// a = the address of a quadword (integer pointer) 
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// b the address of a quadword (integer pointer) 

// c the address of a quadword (integer pointer) 

// 

// Result: c = b - a; 

/////////////////////////////////////////////////////////////////// 
void routine subm (int, int *a, int *b, int *c) 
{ *(int64 *) c = *(int64 *) b - *(int64 *) a;} 

/////////////////////////////////////////////////////////////////// 



IS NEGATIVE ISA REAL NEGATIVE 



/////////////////////////////////////////////////////////////////// 
// Brief description: test the sign bit of a real 
// Expected: a = pointer to a real 
// Result: true if negative 

/////////////////////////////////////////////////////////// / i / / / / / / 
inline int routine is__negative_real (real * a) 
{ if (a->ip) 

return a->ip < 0; 
else 

return a->fp < 0;} 

/////////////////////////////////////////////////////////////////// 



IS ^NEGATIVE IS A 64-BIT INTEGER ARRAY NEGATIVE 



/////////////////////////////////////////////////////////////////// 
// Brief description: test the sign bit of a 64 -bit integer array 
// Expected: a = pointer to start of a 64 -bit integer array 
// (or a real) . 

// n length of array (number of 64 -bit integers) 

// Result: true if negative 

/////////////////////////////////////////////////////////////////// 
inline int routine is_negative (void *a, int n) 

{ return (int) ( { ( (uint64 *) a) [n - 1] » 63) & 1); /* shift the 64-bit integer down*/} 

/////////////////////////////////////////////////////////////////// 



n IS A REAL LESS THAN ANOTHER 



'S5» / // // // // // // // // // // // // // // // // // // // // // // // // // // // // // // // // // 

f*| // Brief description: is a < b 

// Expected: a = pointer to a real 

^ // b pointer to a real 

■fl If Result: result a < b 

^ / // // // // // // // // // // // // // // // // // // // // // // // // // // // // // // // // // 

inline int Iss (real * a, real * b) 

f*sj { if (a->ip < b->ip) 
rZ return 1; 

Ul if (a->ip > b->ip) 
g return 0; 

,,-cw. return a->fp < b->fp;} 

LJ . 1 ( 1 f i f 1 * f f i ' * i ' t i i t i i t i * i t i i i i i i t i i i i t i i i i i t t i i t i i f > f i f i / t i i i t I i i t i 



u NEGATE A 64-BIT INTEGER ARRAY 

/////////////////////////////////////////////////////////////////// 
M: // Brief description: negate a 64 -bit integer array (make negative positive or 
fl! // positive negative) 

Lj // Expected: a = pointer to start of a 64-bit integer array 
laJ // (or a real) . 

ftl II n length of array (number of 64 -bit integers) 

// Result: void 

/////////////////////////////////////////////////////////////////// 
inline void routine negate (void *a, int n) 
{ int i; 

((int64 *) a) [0] = - ( (int64 *) a) [0] ; /* two's compliment*/ 
for (i = l; i < n; i++) 

((int64 *) a) [i] = ~((int64 *) a) [i] ; /* one's complement*/} 

t ' / L L t L L L I * ' I L L * L L I L L f t t t t / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / 



COPY A 64-BIT INTEGER ARRAY 



/////////////////////////////////////////////////////>///////////// 

// Brief description: inline copy a number of 64 -bit integers from 'b' to 'a' 

// (the pointers can point at an integer array or a real) 

// Expected: a = pointer to start of a 64 -bit integer array 

// (or a real) . 

// b pointer to start of a 64 -bit integer array 

// (or a real) . 

// n length of array (number of 64 -bit integers) 

// Result: void 

'////////////////////////////////////////////////////////////////// 
inline void routine cpy (void *a, void *b, int n) 
{ int i; 

for (i = 0; i < n; i++) 

{(int64 *) a) [ij = ((int64 *) b)[i];} 

I ' f ' I i > I f I I J / I I ' I t i I I t I I i f I i f I t t I t I I t 1 I I } l f i l f t l I I I j ! I / I t j l t ! I I / t I ! 



h SHIFTUP BY ONE BIT A 64-BIT INTEGER ARRAY 

///////////////////////////// ;//////////////////////////////////// / 
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// Brief description: shiftup by one bit a 64-bit integers array 
// (the pointers can point at an integer array or a real) 
// Expected: a = pointer to start of a 64 -bit integer array 
// (or a real) . 

// n length of array (number of 64 -bit integers) 

// Result: void 

/////////////////////////////////////////////////////////////////// 
inline void routine shiftup (void *a, int n) 
{ int c, carry = 0; 
int i; 

for (i = 0; i < n; i++) 

{ c = (int) (({uint64 *) a) [i] » 63) & 1; /* shift the 64-bit integer down*/ 

( (int64 *) a) [i] «= 1; 
( (int64 *) a) [i] |= carry; 
carry = c; 

iff )\ / // // // // // // // // // // // // // // // // // // // // // // // // // // // // // // / 



// SHIFTDOWN BY ONE BIT A 64-BIT INTEGER ARRAY 

/////////////////////////////////////////////////////////////////// 

// Brief description: shiftdown by one bit a 64 -bit integers array 

// (the pointers can point at an integer array or a real) 

// Expected: a = pointer to start of a 64 -bit integer array 

// (or a real) . 

// n = length of array (number of 64 -bit integers) 

// Result: void 

/////////////////////////////////////////////////////////////////// 
inline void routine shiftdown (void *a, int n) 
{ int64 c, carry = 0; 
int i; 

for (i = n - 1; i >= 0; --i) 

{ c = (int) (({uint64 *) a) [i] ) & 1; 

( (uint64 *) a) [i] >>= 1; /* unsigned shift down*/ 
( (uint64 *) a) [i] | = (carry « 63); 
carry = c; 

}} 



LSS COMPARE FUNCTIONS FOR AN ARRAY OF 64-BIT INTEGERS 

'/////////////////////////////////////////////////////////////////// 

fl // LSS LESS THAN COMPARE AN ARRAY OF 64 -BIT INTEGERS 

|H // Expected: 

// a address of 64 -bit integer array 

* II b address of 64 -bit integer array 

O // n number of 64 -bit integers 

Ui II Result: TRUE if (a < b) else FALSE; 

w*? / // // // // // // // // // // // // // // // // // // // // // // // // // // // // // // // // // 
y» inline int routine lss (void *a, void *b, int n) 
ssk 2 { int i; 

ijj if (({int64 *) a) [n - 1} 1= ( (int 64 *) b) [n - 1]) 

Lj return ( (int64 *) a) [n - 1] < ((int64 *) b) [n - 1] ; /* signed compare*/ 

f§s for (i = n - 1; i >= 0; --i) 

{ if (<(int64 *) a) [i] i= ( (int64 *) b) [i] ) 

return { (uint64 *) a) [ij < ( (uint64 *) b) [i] ; /* unsigned compare*/ } 

return 0; /* false equal*/} 

/////////////////////////////////////////////////////////////////// 



u LEQ LESS THAN OR EQUAL COMPARE TWO REALS 

/////////////////////////////////////////////////////////////////// 
// Expected: 

// a address of real 

// b address of real 

// Result: TRUE if (a <= b) else FALSE ; 

/////////////////////////////////////////////////////////////////// 
inline int routine leq (real * a, real * b) 
{ if (a->ip <= b->ip) 
return 1; 
if (a->ip > b->ip) 

return 0; 
return a->fp <= b->fp;} 

/////////////////////////////////////////////////////////////////// 



// LEQ LESS THAN OR EQUAL COMPARE AN ARRAY OF 64-BIT INTEGERS 

/////////////////////////////////////////////////////////////////// 
// Expected: 

// a address of 64 -bit integer array 

// b address of 64 -bit integer array 

// n number of 64 -bit integers 

// Result: TRUE if (a <= b) else FALSE; 

/////////////////////////////////////////////////////////////////// 
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inline int routine leq (void *a, void *b, int n) 
{ int i; 

if (<(int64 *) a) [n - 1] != ( (int64 *) b) [n - 1]) 

return ((int64 *) a) [n - 1] < ( {int64 *) b) [n - 1] ; /* signed compare*/ 
for (i = n - 1; i >= 0; --i) 

{ if (((int64 *) a) [i] 1= ((int64 *) b) [i] } 

return ( (uint64 *) a) [i] < ( (uint64 *) b) [i] ; /* unsigned compare*/ } 

return 1; /* true equal*/} 

//////////////////////////////////////// / / / // // // // // / / / iii i i ii i i i i 

„ GEQ GREATER THAN OR EQUAL COMPARE AN ARRAY OF 64-BIT INT EGERS 

I I l l i I I i f I t I I I I t I l I l I I I I I l I I / l I I i 1 I I I I i I I I I / t l I I I ! I I I I I I I i ! I < I I I I I I 

If Expected: 

// a address of 64 -bit integer array 

// b address of 64 -bit integer array 

// n number of 64 -bit integers 

// Result: TRUE if (a >= b) else FALSE; 

/////////////////////////////////////////////////////////////////// 
inline int routine geq {void *a, void *b, int n) 
{ int i; 

if ({{int64 *) a) [n - 1] != ( (int64 *) b) [n - 1]) 

return ((int64 *) a) [n - 1] > ((int64 *) b) [n - 1]; /* signed compare*/ 
for (i = n - 1; i >= 0; --i) 

{ if {{(int64 *) a) [i] != ({int64 *) b) [i] ) 

return ( (uint64 *) a) [i] > ( (uint64 *) b) [i] ; /* unsigned compare*/ } 

return 1; /* true equal*/} 

i i i i i i i t i i i i i i i i i i i i i i i i i t i i i t i i i i i i i t / i t i i i i i i t i i i i t i f i i i i i i i i i i i i 

b SET TO ZERO A 64-BIT INTEGER ARRAY _____ 

///////////// i //////// i i i i ii i i iii i iii i ii t i i i i i i i i i i i i i i i i i i i i i i i i i i 

II Brief description: set to zero a 64-bit integer array 
// (the pointers can point at an integer array or a real) 

// Expected: a = pointer to start of a 64 -bit integer array 
// (or a real) . 

// n length of array (number of 64 -bit integers) 

// Result: void 

i / i i i i i i i i i / i i / i i i i i i i i i i i i i i i i i i i i i i i i i i / i i i i i / i i i i i i i i i i i i / i i / i / i 
inline void routine set_zero (void *a, int n) 
{ int i; 

for (i = 0; i < n; i++) 
( (int64 *) a) [i] = 0; } 

//// / /////////////////////// i // i i i i i i i i i i / i / i i i i i i i i i i i i i i i i i i i i i / i 

// TEST IF A 64-BIT INTEGER ARRAY IS ZERO 

i //// i ///// i i i i i / i i / i // i i i i i / i i i i i i i i i i i i i i i i i / i i i i i i i i i i i i i i i i i i i i 

II Brief description: test if a 64 -bit integer array is zero 

// (the pointers can point at an integer array or a real) 

// Expected: a = pointer to start of a 64 -bit integer array 

// (or a real) . 

// n length of array (number of 64 -bit integers) 

// Result: void 

i i i i i i i i i i i i i i i i i i i i i i i i i i i i t i i i i i i i i i i i / i i i i i i i i i i i i i i i i t i i i i i i i i i 
inline int routine is_zero (void *a, int n) 
{ int i; 

for (i = 0; i < n; i++) 
if ( ( (int64 *) a) [i] ) 
return 0; 

return 1; /* yes zero*/} 

I I I I I i f I I f I f i l t f i t f t i t l I f I I I f I I I I I i t I I I I I I t t I l I I I I I I I I l I I l l I ! I l I I I i 

h ADD64 ADD 2 64-BIT INTEGERS AND PLACE RESULT IN A THIRD 

i ////// i ////// i / / i / i i / i i i i i // i i i i i i i i i i i i / i i i i i i i i i i i i i i i i i i i i i i i i i 

II Brief description: add 2 reals and place result in a third 
// Expected: 

// a = the address of a 64 -bit real 

// b the address of a 64 -bit real 

// c = the address of a 64 -bit real 

// 

// Result: c = a + b; 

i i i i i i i i / i i i i i i i i i i i i i i / / i i i i i i i i t i i i t i i i i i i t i i i i i i i i t i i i i / i i i i i i i / 
void routine add64 (real * a, real * b, real * c) 
{ *(int64 *) c = *(int64 *) a + *(int64 *) b; } 

i i i i i i i i i i t i i i i i i i i i i i i i / i i i i i i i i i / i i i i / i i i i i i i i i i i i i i i / i i i i i / i i i i i 

n SUB64 SUBTRACT 2 64-BIT INTEGERS AND PLACE RESULT IN A THIRD 

i i / t i i i i i i i i i i / i i i i / i i i i / i / i i / i / i / i i i i i / i / i i i i i i i i i / i / i i i / i i i i / i i i i 
II Brief description: subtract 2 reals and place result in a third 
// Expected: 

// a the address of a 64 -bit integer (void pointer) 

// b = the address of a 64 -bit integer (void pointer) 

// c = the address of a 64 -bit integer (void pointer) 

// 
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/ / f / / / ///////////////////////////////// ' ' ' ' ' ' ' 11 ' ' ' / ' / ' 1 ' 1 1 1 ' / 7 7 ; 
void routine sub64 (void *a, void *b, void *c) 

{ *(int64 *) c - *(int64 *) a - * (int64 *) b;} ,>,,,,,>,,,,,, i t t 

1 tt ADD A 64-BIT INTEGER ARRAY (LEN GTH OF TWO) 

1 i i t , i i , t } i f t t i i i t i i t i i i t t t t i i i t i i i i i i i i i i i i i > * ' i i 1 1 1 1 1 1 ' ' ' 1 ' ' ' 1 ' 1 ' 1 
If Brief description: inline add a length of two 64 -bit integers a = a + b 
// slightly faster for a 'normal' add of a 128-bit integer 
// (the pointers can point at an integer array or a real) 
// Expected: a = pointer to start of a 64 -bit integer array 
// (or a real) . 

// b pointer to start of a 64-bit integer array 

// (or a real) . 

w / / R ? S / U / V / / / V / OX / d / ////////////////////////////////'/ /////// / / / / / / ' / ' / ' / 
inline void routine add2 (void *a, void *b) 
{ uint64 was; 

int64 carry = 0; 

was = ( (uint64 *) a) [0] ? 

((uint64 *) a) [0] += { (uint64 *) b) [0] ; 

carry = ( ( (uint64 *) a) [0] < was) j 

{(uint64 *) a) [1] += ( (uint64 *) b) [1] + carry;} ,,,,,,,,,/,/,/,,/ y 

//////////////////////////////////////////////// 'I L 1 i 1 i 1 1 1 1 1 1 1 1 1 ' 1 1 

1 h SUBTRACT A 64-BIT INTEGER ARRAY (LE NGTH OF TWO) 

X l I { t t l ! I l I I l f f I ! I t i I i I I I I I I I I I l I I I I l I I t I I I I i I I i i I ! I * f I 1 1 I I 1 1 1 1 1 1 1 ' 1 

If Brief description: inline subtract a length of two 64 -bit integers a = a - b 

// slightly faster for a 'normal 1 subtract of a 128-bit integer 

// (the pointers can point at an integer array or a real) 
\& II Expected: a = pointer to start of a 64 -bit integer array 
isa. // (or a real) . 

y // b pointer to start of a 64 -bit integer array 

C3 // (or a real) . 

• // Result: void 

t i t t t i t i i t t t t i t t t i i t i t t t t t t i t t t t t t t i t i t t i i t t t t t i t t / i t t t i i t t t i § t i i i i 
yl inline void routine sub2 (void *a, void *b) 

{ uint64 was; 
^ int64 carry = 0? 

D was = ({uint64 *) a) [0] ; 

( (uint64 *) a) [0] -» ( (uint64 *) b) [0] ; 
ws carry = ( ( (uint64 *) a) [0] > was) ; 

s ((uint64 *) a) [1] -= ( (uint64 *) b) [1] + carry;} ,,,,,,,, y , / > 

^ ///////////////////////////////// / / / / / / / / / / / / / / / / / / / / / / / / t i i i i ' 1 i 1 ' 



n SCALEUP A 128-BIT INTEGER BY SCALE . 

i i , i , i i i i i i i t i i t i i t i i t t i i i i i i i f i i i i i i i i i t t i i i i i 7 / / / / / i i / / / / 7 / i i i i i i 
II Brief description: Multiply the 128 -bit integer by SCALE 
// Expected: 

// _a = address of 64 -bit integer array (_a * Jo) 

// Result: void 

/ z / / / / / / / / / / / / / / / / / / / / / t i / / / / / / / / / / / / / / / / / / / / t / / / / i i i i i i i 1 i i i 1 i i i { 1 
inline void routine scaleup (void *_a) 
{ int64 a[] = { 0 r 0 }; 

int64 b[] = { SCALE, 0 }; 

int64 c[] = { 0, 0 }; 

int na = 0; 

int i; 

cpy (a, _a, 2) ; 

if (is_negative (a, 2)) 

{ /* negative*/ 

negate (a, 2); /* make positive*/ 

na - 1; } 
if (lss (a, b, 2)) 

{ /* swap*/ 

int 64 temp [2] ; 
cpy (temp, b, 2) ; 
cpy (b, a, 2) ; 
cpy (a, temp, 2) ; } 
for (i - 0; i < 64 && (b[0] || b[l]); i++) 
{ if (b[0] & 1) 

add2 (c, a) ; 
shiftdown {b, 2) ; 
shiftup (a, 2) ; } 
if (na == 1) 

{ negate (c, 2); /* make positive*/ } 

cpy ( a . c , 2 ) ; } 

,/////-////////////////////////////////////////////////''/'''''''''' 
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0 SUBTRACT A 64-BIT INTEGER ARRAY . 

////////////////////////////////////////// /// ' ' ' ' 11 / ; / / ' 1 < 1 1 1 1 1 ' ' 1 1 

fl Brief description: inline subtract a number of 64 -bit integers a = a - b 

// (the pointers can point at an integer array or a real) 

// Expected: a = pointer to start of a 64-bit integer array 

// (or a real) . 

// b pointer to start of a 64-bit integer array 

// (or a real) . 

// n length of array (number of 64 -bit integers) 

W / W^lV / / / V / 01 / d / I I 1 I I t i I I I t I ! I I t I I I t I I I I t I I t I I I I I I I < I I t I I I I ' I ' 1 I 1 ' 1 1 1 

inline void routine sub (void *a, void *b, int n) 
{ uint64 was; 

int 64 carry = 0; 

int i; 

for (i - 0; i < n; i++) 

{ was = ( (uint64 *) a) [i] ? 

((uint64 *) a) [ij -= ( (uint64 *) b) [i] + carry; 
carry = ( ( (uint64 *) a) [i] > was) ; 

t ,i )) /, / /// f // / / t / t /// / / / t / // t ////// t t //'/'' / ' ii LJ. i L I' L i 1 ! 1 1 1 1 1 1 1 1 1 

a ADD A 64-BIT INTEGER ARRAY „ 

////////,///////////// i // t // i i i / / i // ' i ' ' ' i ' ' i i ii i 1 / ' 1 { * ' / ' f ' ' 1 1 1 1 ' ' 

II Brief description: inline add a number of 64-bit integers a - a + b 

// (the pointers can point at an integer array or a real) 

// Expected: a = pointer to start of a 64 -bit integer array 

// (or a real) . 

// b pointer to start of a 64-bit integer array 

// (or a real) . 

// n length of array (number of 64 -bit integers) 

i / R ? S / U / V / / 71*?' 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 { 1 1 1 1 > 1 1 1 1 1 1 1 

inline void routine add (void *a, void *b, int n) 
{ uint64 was; 

int 6 4 carry = 0; 

int i; 

for (i = 0; i < n; i++) 

{ was = { (uint64 *) a) [i] ; 

( (uint64 *) a) [i] += ( (uint64 *) b) [i] + carry; 
carry = ( { (uint64 *) a) [i] < was) ; 

^ I l I I I I t I I I I I I 1 I I t I ! I I i ! It f I t I ! I I I I I I I / I I I I I I I I 1 I I I I I 



I I I I I t I I t I I I I I I 



u SCALEDOWN A 128-BIT INTEGER B Y SCALE 

t I I I I { I t I I l I I I f I I l t I t i I I I I I I I I I I J l I I ! I I I / I I II I f t I I t I I f I t I 1 ! I I I I I I I I 

II Expected: 

// _a = address of 64 -bit integer array 

// Result: void 

/////////////////////////////// / / // // / / / tit i i i i i //////// i i i i i i i i i i i 
inline void routine scaledown (void *_a, int n) 
{ uint64 a[] = ( 0 ( 0, 0, 0 ); 
uint64 r[] = { 0, 0 }; 

uint64 scale [] = { SCALE, 0, 0, 0 }; /* SCALE must be less than an int 64*/ 
uint64 b[] = { SCALE, 0, 0, 0 }; 
uint64 pow[] = { 1, 0 }; 
int na = 0; 
int i; 

cpy (a, _a, n) ; 
if (!a[lj && Ia[2] && !a[3]) 
{ a[0] /= SCALE; 

cpy (_a, a, n) ? 
return ; } 
if (is_negative (a, n) ) 

{ /* negative*/ 

negate (a, n) ; 
na = 1; } 
/* shift »b' until it is greater than 'a'*/ 
for (i = 0; i < 128 && leq (b, a, n) ; i += 2) 
{ shiftup (b, n) ; 

shiftup (pow, 2) ; 
shiftup (b, n) ; 
shiftup (pow, 2) ; } 
if (is_zero (pow, 2)) 

{ runtime$error (MSG_FL0AT0VF) ; /* overflow*/ } 

while (leq (scale, a, n) ) 
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{ /* while still a divisor*/ 

while {lss (a, b, n) ) 
{ shiftdown (b, n) ; 

shift down {pow, 2) ; } 
if (!a[l] && ia[2] && !a[3] ) 
{ /* down to a 64-bit integer?*/ 

a[0] /= SCALE; 
add2 (r, a) ,- 
break; } 
sub (a, b, n) ; /* subtract largest power*/ 

add2 (r, pow) ; /* increment result*/ } 

if (na == 1) 

{ negate (r, 2) j } 

cpy (_a, r, 2) ; } 



/////////////// 



I I I I I ( I I I I I ! I I t I I t i i I I I I I I 1 I I I I I 1 I I I I I I If ! I I i I' I I I ' ' 



REALIP GET INTEGRAL PORTION OF A REAL 

I t I [ I I I I I t I I I I I I l I t I I t t I I I I I t I I I I I I I I i I I ! I I I i I I I I I * I I I I l I I l I * 1 1 ' 1 1 1 

ff Brief description: get the integral portion of a real removing the fractional part. 
// Expected: 

// in = address of input real 

// out = the address of a 64-bit real {void pointer) 

// Result: void : out = integral portion of in 

/ / // // // // // // // // // // // // // // // // // / / / // // // // / / / / / / / / / / / / / / / / / / / / 

void routine realip (real * in, real * out) 

{ out->fp =0; /* no fractional part*/ 

out->ip = in->ip; /* just the integral part*/} 

///////////////////// / / ////// / ////// / / / / / / / / / / / / / / / / / / / / / / i i i i i i i i i 



WHOLE INT GET INTEGRAL PORTION OF A REAL 



///////////////////// /////////// / / // // // / // // // / // / / / / / ! I I I I l I I t I l I 

p5 // Brief description: get the integral portion of a real removing the fractional part. 
2; // Expected: 

O // in = address of input real 

,fi // out = the address of a 64 -bit real {void pointer) 

~"Z ff Result: void : out = integral portion of in 

yl ////////////////////////////////////// i /// i // i / i i i i i i i i i i i i i i i i i i i i 

void routine whole_int (real *in, real *out) 
-aw { realip (in, out);} 

I I j t I I [ I I f I I f I t f t I I I I t I I I f I I I I I I ! ! I I t I I I I I t I ( I I t I * I I t * I ( I I * 1 1 < ! 1 1 1 f 



h CVTLR CONVERT LONG TO A REAL 

y?i ////////////////////////// i /// i /// / i ii i ii i ii i i i i i i i i i i i i i i i i i i i i i i i 

If Brief description: convert an integer to a 64 -bit real 
// Expected: 

\g£ ff a = input integer 

yj // b = the address of a 64 -bit real 

s"T // Result: 1 : b = a 

jp& /////////////////////////////////////////////////////////////////// 

f*-$ int routine cvtlr (int _a, real * b) 

' : jf { real a = _a; 
Q *b = a; 
ff} return 1 ; } 

5 ^ /////////////////////////////////////I///////// I / I ' I II I I II I / / 'It I I I 



n CONVERT REAL TO A LONG 

////////////////////// ///////// / ////////// / / J / / / / / / / / / 7 / / i i i i i i i i i i 

If Brief description: convert a real to a long with rounding 
// Expected: 

// a address of input 64 -bit real 

// b = address of output integer 

// Result: 1 : a = rounded real b 

t i i i i i / t i / t t i t t / i f i i t i i i i i i t t t i t i i i i i i i i i t i i i * ( i i ' i i i i i 1 * i t * i i * 1 * * i 
int routine cvtrl (real *r, int *i) 

{ *i = (int) r->ip; // get integer part 

if (r->fp i= 0 // Non-ZERO fraction? 

&& (r->fp >= SCALE / 2 // and we have to round it 
| | r->fp <= - (SCALE/2) ) ) 

{ if (*i < 0) // round down if negative 

*i = *i - 1; 
else 

*i=s*i + l; // round UP if zero or positive 

} 

return 1 ; } 

/ / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / i i i i i i i 



GETINT GET AN INTEGER FROM A REAL 



/ ///////////// / ///////// / ////// / ////// / / / / / / / / / / / / / t / / 7 7 / 7 / / / i i i i i i 
ff Brief description: Get an integer from a real 
// Expected: 

// n reference to a real 

// Result: integer value 
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/////////////////////////// i i i f i i i i i i i i i i i i i t i i i i i i i i t i t i t 1 1 f i i 1 f i t 
int routine getint (real & n) 
{ int i; 

cvtrl (&n, &i) ; 

return i ; } 

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*%%*%%%*** 
// I P 

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
// Brief description: Truncates a real number at the decimal point and 
// returns the integer portion. 
// Expected: in - passed real number 

// Result: out - the integer portion of a real number 
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%^ %%%%%%%%% 
void routine ip (int *out, real * in) 

( *out = (int) in->ip; /* trunctate a int 64 down to an int*/} // ip 

/ / / / / / / / / / / / / , ' / i / / / / / / / / / / / / j i / / / / / / / / / / / / / / / / / ' / / i 1 1 1 i 1 1 ' 1 1 1 1 ' 1 1 s 1 

n CMP64 COMPARE TWO 64-BIT INTEGERS 

/ / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / i * t ' 1 1 1 1 1 ' 1 1 1 1 1 { 1 
II Brief description: compare two 64 -bit integers 
// Expected: 

II a = address of 64 -bit integer 

// b = address of 64 -bit integer 

// Result: 0 = equal 

// -1 = less than 

// 1 = equal 

///////////////////////////////////////////////////'/''/' ' ' ' 1 1 f f 1 ' 1 
int routine cmp64 {void *a, void *b) 
{ if (*(int64 *) a *(int64 *) b) 

return 0; 
if (*(int64 *) a < * (int64 *) b) 

return -1; 

/ / ^lj^±Llj_ i i i i i I i i i i i i i t i i I i i i i i i i i i l i i i i 1 i i 1 1 1 > 1 ' i 1 ' 1 1 1 ' 1 ' 1 ' 1 1 1 1 1 1 ' 



„ DHEX64 DUMP A 64-BIT INTE GER IN HEX - 

/////////////////// ///////////////////////////'// i i it i 1 1 f 1 1 1 1 t 1 1 1 ' 1 
II Brief description: Dump a 64 -bit integer in hex. This is a debug tool 
// for 64-bit arithmetic. 
// Expected: 

// buf - dump hex characters to this buffer 

// _s = start address of 64 -bit quadword 

// Result: void 

; / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / ' / ' / / ' / ' t ' ' ' * 1 

static void routine dhex64 (char *buf , void *_s) 
{ unsigned int *s = (unsigned int *) _s; 

sprintf (buf, "%08X%08X» / s [1] # s[0]);} ,,,,,,,, 

////////////////*//////////////////// / It! t / / I I I t I f I I I I I I I I I I I I I t f I I I 



n DHEX96 DUMP A 96-BIT INTEGER IN HEX 

I , , I l l , , f I I t I } I I I I I l ! ! I I I I t t I l I I I I ! I l I I I I i I l t I I i I I I i I l I I I * > * } 1 1 1 1 1 1 

II Brief description: Dump a 96 -bit integer in hex. This is a debug tool 
// for 96-bit arithmetic. 

// Expected: 

// buf = dump hex characters to this buffer 

// _s = start address of 96 -bit quadword 

// Result: void 

f/fft//if//f/f///fff(fftt/i/f/t/ffff/ftffttf( / // / / / / // // // // // // // / 
void routine dhex96 (char *buf , void *_s) 
{ unsigned int *s = (unsigned int *) _s; 
dhex64 (buf , s + 1) ; 

sprintf (buf + 16, "%08X", s[0]);} 

/////////////////////////////////////////// LJ. f i { i 1 / f 1 { / 1 1 ! ' 1 ' ' ' 1 1 1 1 



n DHEX128 DUMP A 12 8-BIT INTEGER IN HEX | 

/////////////////////// //////////////////'/ / / // til I I I J I I I ! I I I I I i I I I l 

// Brief description: Dump a 12 8 -bit integer in hex. This is a debug tool 
// for 128-bit arithmetic. 
// Expected: 

// buf = dump hex characters to this buffer 

// _s = start address of 128-bit quadword 

// Result: void 

/ / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / ' ' ' / ' ' ' ' 
void routine dhexl28 (char *buf, void *_s) 
{ unsigned int *s = (unsigned int *) _s; 

dhex64 (buf, s + 2) ; 

dhex64 (buf + 16, s) ; } 

I l I I f I f I ! I t f l I t I i f I I ! I t i I l t l I t t l I t J f ! I I ! I l ! t i ! t f l I I ! I 1 i t ! 1 I i I t i I t I t 

RTOD CONVERT REAL TO DOUBLE | 



/ / / ? / 7 / / / / / / / / / / / / / / / / / / / / / / ? / / t i / / ( ) / / / / / / / / / f i f i t i i t * i i ' ' 1 i i i 1 < 1 i 
II Brief description: Convert a real to a double. 
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// Expected: 

// a address of real 

// b address of double 

// t Result : t void f ^ f i ^ f f t f t / f / / / / f t f f f f f f f f f t f f f f t ,,,,,,,,,,,,,,,,,,,,/ / 

void routine rtod (real * a, double *b) 
{ double d - a->fp; 

d /= SCALE; 

d += a->ip; 

/ / * b ; 7 / d / \ / / / / / / / / / / / t t / / / / / L_L / / / / L ' / > L > i ' L > ' L ' > L ' ' L I ' 1 > 1 ' ' 1 1 1 > / ' ' 1 11 

• GETINT GET A DOUB LE FROM A REAL ,,,,,,,,,, 

////////////////////////////////////// / rn / / / / / / / / 7 / / / ' * ' ' * > * f 1 1 1 ' ' ' 
II Brief description: Get a double from a real 
// Expected: 

// n address of a real 

, / f^rrrYr, ,,,,,,, , ,,,,,, , ,,,,,,,,,,,,•,,*,>•">>> > < 

double routine getdouble (real * n) 
{ double d; 
rtod (n, &d) ; 

/ / f^y? 1 /^! , 1 1 1 1 1 1 1 1 , 1 1 f , 1 1 1 1 1 * 1 1 * 1 1 1 1 1 1 1 1 1 * 1 1 1 * i t* / ul > < i 1 1 1 1 1 1 1 1 > 1 1 1 

„ DIV128 DIVIDE TWO 128-BIT INTEGERS . 7 __ 

/////////////////////////////////// ////// t * 1 f 1 1 / / , / / / / 1 ' 1 1 1 1 1 ' 1 ' 1 1 1 



II Brief description: Divide two 128-bit integers. This is used for intermediate 
// arithmetic. 
// Expected: 

// _dvend = address of real dividend 

// _dvsor = address of real divisor 

// quote = address of real quotent (division result) 

// rem = address of real remainder 

// Result: dividend/ divisor = quotent, remainder .,,,,,,,,,,,,,,,//////// 
void divl28 (void *_dvend, void *_dvsor, void *_quote, void *rem) 
{ uint64 dvend [] = { 0, 0 }; 

uint64 sdvsor[] = { 0, 0 }; 

uint64 dvsor [] = { 0, 0 }; 

uint64 quoted = { 0, 0 }; 

uint64 pow[] = { 1, 0 }; 

int i ; 

cpy (dvend, _dvend, 2); /* get local copy of dvend*/ 
cpy (dvsor, _dvsor, 2); /* get local copy of dvsor*/ 
cpy (sdvsor, _dvsor, 2) ; /* get local copy of dvsor*/ 
if (is_zero (dvsor, 2)) 

{ /* no divisor ?*/ 

runtime$error (MSG_DIVBY0) ; /* divide by zero*/ } 
/* * shift divisor until it is >= to dividend 

* that is, whilst the divisor is less than the dividend.*/ 
for (i = 0; i < 128 && lss (dvsor, dvend, 2); i += 2) 

{ shiftup (dvsor, 2) ; 

shiftup (pow, 2) ; 

shiftup (dvsor, 2) ; /* shiftup twice to speed up*/ 

shiftup (pow, 2) ; } 
if (is_zero (pow, 2)) 

{ runtime$error (MSG_FLOATOVF) ? /* overflow*/ } 

/* * Here we subtract the largest possible divisor value 

* from the divident each time adding power to the quotent. 

* We are trying to reduce the divident as quickly as we can.*/ 
while (leq (sdvsor, dvend, 2)) 

{ /* while still a divisor*/ 

while (lss (dvend, dvsor, 2)) 
{ /* follow the divident down in value*/ 

shift down (dvsor, 2) ; 
shiftdown (pow, 2) ; } 
if (! sdvsor [1] && ! dvend [1]) 
{ /* down to just 64-bit integers ?*/ 

uint64 n[] = { 0, 0 }; 
n[0] = dvend[0] / sdvsor [0]; 

dvend[0] -= n[0] * sdvsor [0]; /* remainder*/ 

add2 (quote, n) ; 

break ; } 
sub2 (dvend, dvsor) ? /* subtract largest dvsor*/ 
add2 (quote, pow) ; /* increment quotent*/ } 
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if (rem) 

cpy (rem, dvend, 2) ; /* remainder*/ 
cp y/ (.quote ; quote, 2)^ , /^return result;/ } f ,,,,,,,,,,,,,,,,,,,,,,,,,,, ,,,,,,, 

I „ MUL128 MULTI PLY TWO 128-BIT INTEGERS ,,,,,,,,,,,,,,> 

1 /////////////////////////////////// ' < > ' if ' 1 f * > J 1 } 1 '! ! 1 ' ' ' ' 1 ' ' 1 ' ' ' 1 ' 
II Brief description: Multiply two 128-bit integers. This is used for intermediate 

// arithmetic. 

// Expected: 

// _a = address of 64-bit integer array (_a * _b) 

// ~b = address of 64-bit integer array (_a * _b) 

// r address of 64 -bit integer array (result) 

// ^ Result : t void ^ f t f f t f ( t f f f f ( f f f f ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,/,// / 
void mull28 (void *_a, void *_b, void *_c) 
{ uint64 a[] = ( 0, 0 ( 0, 0}; 

uint64 b[] = { 0, 0, 0, 0}; 

uint64 c[] = ( 0 ( 0, 0, 0}; 

int i; 

cpy (a, _a, 2) ; 
cpy (b, _b, 2) ; 

/* * make 'b 1 the smaller of the two*/ 
if (Iss (a, b, 4)) 

{ int 64 temp [4] ; 

cpy (temp, a, 4) ; 

cpy (a, b, 4) ; 

~ for (i = 0; i < 128 && ! is_zero (b, 4); i++) 

p if (b[03 & 1) 

add (c, a, 4) ; 
shiftdown (b, 4); /* --> shift down*/ 

d3 shiftup (a, 4); /* <-- shift up*/ } 

/ / J P 7 , full ill ±ll i t t t i t i i t t i t t i t t t t i * t it i_ L (-L / / L-L L i f t L ' f 1 ' ' i 1 1 ' 1 1 ' 1 ' 1 ' 1 ' 
O I n DIVR DIVIDE T WO 128-BIT REALS . _— 

:Z X ( I , , t f f , I l I I I t I { I 1 I t t ! t ( t l I t t I I I 1 ( I l 1 l t { i ! I 1 ! I t t ! I I I i I ' I I t i l t 1 f < I t > 

US // Brief description: Divide two 128-bit integers. This is used for intermediate 
rs // arithmetic. 
*^ // Expected: 

^ // _dvend = address of real dividend 

yj // _dvsor = address of real divisor 

lH If quote = address of real quotent (division result) 

!Z // rem = address of real remainder 

IU II Result: dividend/divisor = quotent, remainder t , t * t t 

y void divr (real * dvend, real * dvsor, real * quote, real * rem) 
fly { real q; 

r96 divide (dvend, dvsor, &q) ; 

if (rem) *rem = 0; 

/ / Vj^L-LJ. V) f > i f i t ' i > > i > < ' > ' i ' > ' t i > > ' * i i f 1 1 ' > > 1 1 1 < ! ' 1 ' 1 1 1 1 1 1 1 1 1 ' 1 * * * f 



0 REAL::OPERATOR-() NEGATE A REAL 

z / ////////////////////////////////////////// ///// ' ' / / ! 1 1 ' ' 1 ' ' ' 1 1 ' ' 1 
II Brief description: negate the current real as in -real 
// Expected: The current real class 

i 1 , , fTfY'i i / z*?* 1 / ////////////////////////////////'//' ' ' 1 ' 11 1 1 1 1 f { 1 ' 1 1 ' 
real real :: operator - () 
{ real r - *this; 
negate_real (&r) ; 

/ / f e / t y? l / r / / ///////////////////////////// / ' LS-1 LJ. L <-L LU f—L * L ' L L t i 1 1 1 1 1 / t > 

/, REAL::OPERATO R+=() ADD A REAL , , , , , 

1 f j 7 , f i , , i , / i i , i i i i i i i i i i i i i i i / / / / / / / / / / / / / / / / / / i i < i ! ' ' 1 1 ' 1 1 1 { 1 1 ' 1 1 
II Brief description: add a real to ourselves as in real += real 

// Expected: The current real class 

// ^ Result : f t real^ reference 7 , /7 , ; ; y , ; ;/ /7 //7 7// ; /7 ;/ 7/ 7 y/ //; //; 7/ ; /,///,, , 
real &real :: operator += (real r) 
{ /* real += real*/ 

addr (this, &r) ; /* add real*/ 

/ / fyy 37 / 1 /*^ 1 / 8 /^/ 1 1 1 1 i 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 f 1 1 1 1 1 1 1 1 * i * * * f 1 1 { 1 * 1 1 1 > > { 1 * * 1 1 * * 1 

» REAL:: OPE RATOR++Q INCREME NT A REAL 

I , , t t ! I , t I 1 I I I t t { t I I t I I l t i I t I I I t I I I t 1 I I ' I I t I' ^ I' 11 1 1 1 1 1 / ' 1 1 / 1 1 1 ' 1 ' 

II Brief description: increment a real by one as in real++ 
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// Expected: The current real class 
// Result: real reference 

I f i 1 ! I I t ! I t I t f i i t I l t t l I t ! I t f I I I I i I I I I I t i l I t I I I J I t i I I I I I I t I * I 1 * 1 ! 1 I * 

real &real : :operator++ (int) 
{ /* real++*/ 

ip++; 

return *this;} 

///////// / tt tttiifiitiifiiiiititittitttiiii i / / f t / / / // // / / // // // // // 

i REAL::OPERATOR-=() SUB TRACT A REAL ^ 

////////////// ////////////////////////////////////// i ii ti i i 1 i 1 > 1 ' i ' 

I I Brief description: subtract a real from ourselves as in real - = real 
// Expected: The current real class 

// Result: real reference 

//////////////////////////////////////// / / / / / / i i i i i i i i 1 i i i i i i i ' i 1 ' 1 

real &real :: operator -= (real r) 

{ /* real -= real*/ 

subr {this, &r) ? /* subtract two reals*/ 

return *this;} 

/ / // // // // // // // // // // / / / / / /////// / / // / i i t i i i i f i t i t i iti i t fi t i iii i i i 

u REAL:: OP E RATO R-Q DECREMENT A REAL 

////////////////// ////// / / / // / / // / ////// i ////// / / / i i i / i / / / / / / / / / // / 
// Brief description: decrement a real by one as in real-- 
// Expected: The current real class 
// Result: real reference 

/////////////////////////////////////////////////////////////////// 
real &real :: operator-- (int) 
{ /* real--*/ 

ip--; 

return *this; } 

//////////////// / / / // // // // / / // / / // / / / // // // / / / / / / / / / / / / / / / / / / / / / / / 
// REAL::RE AL() CONS T RUCT A REAL 

///////////// ///////////////////////////////////// / f i i i i i ii i t t ii i i i 
II Brief description: construct a real and initialize it to zero 
// Expected: The current real class 
// Result: real reference 

///////////////////////////////////////// / / / / / / tit t i i / / / / / // // // // / 
real : : real ( ) 
{ ip = 0; 
fp - 0;} 

i t I I I l I I l i I 1 I I I I I i I I l I l I f I l I t I t t I I I i i t I 1 I I t I I I I I 1 I I I I I ! I I i I I I ' I I f l ' 

u REAL::REAL() C ONSTRUCT A REAL FROM AN INTEGER 

7 7 7 7 7 7 7 7 7 I 7 7 7 7 7 7 7 7 7 1 7 7 7 7 7 7 1 7 7 7 7 7 7 7 7 7 7 T~l 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 

/ / Brief description: construct a real from an integer 
// Expected: The current real class 
// Result: real reference 

/////////////////////////////////// / i i i i i i / / / / / / / / / / / / / / / // // // // // 
real: :real (int i) 

{ /* construct a real from an integer*/ 

ip = i; 
fp = 0;} 

/ / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / i i i i i i i i i i i i i i i i 



n REAL::REA L() CONSTRUCT A REAL FROM A 64-BIT INTEGER 

/////////////////////// i /////// i /// i ////// i // i i i ' i i i i i i ' i ' ' i i i i i ' ' i 
II Brief description: construct a real from a 64-bit integer 
// Expected: The current real class 
// Result: real reference 

//////////////////////////////////////////////////////////////////' 
real:: real (int64 i) 

{ /* construct a real from a 64 -bit integer*/ 

ip = i; 
fp - 0;} 

///////////////////////// i /// i /// / ////// i /// i /////// / i i i i / i ' i / ' ' s i ' 



h REAL::RE AL() CONSTRUCT A REAL FROM A DOUBLE 

"~ 7 7 7 7 7 7 7~"7 7 7 7 7 7 7 7 7 7 7 7~~7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7^ 7 7 7 7 7 7^ 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7^ 7 T~ 
II Brief description: construct a real from a double 
// Expected: The current real class 
// Result: real reference 

/////////////////////////////////////////////////////////////////// 
real:: real (double d) 
{ dtor (Scd, this) ; } 

ffft/ftfitfttftfffi/t/f/fffffftif/ftfffttfttftfftfftftt/ftffttftff/ 
II real (whole, fraction) 

/////////////////////////////////////////////////////////////////'/ 
real:: real (int64 whole, int fract) 
{ ip = whole; 

fp - fract; } 

t t f t i i / t t / i i / i i i { i i i i i t / / ; i t i i i i i f t t / / i i / i i / i i i i t / i t t i t t i t / i i i t i < i i 



n REAL::OPERATOR-() SUBTRACT A REAL FROM ANOTHER REAL 

//////////////// i ////////////////////////////////////////// i i f / / / / / 
// Brief description: subtract a real from another real 
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// Expected: The current real class 
// Result: real reference 

/////////////////////////////////////////////////////////////////// 

real real :: operator - (real i) 

{ /* real - real*/ 

real r = *this; 

return r - = i ; } 

I I 1 I I t I I l ! i I i i I t ! I I ! I t I ! I I I I t I I I l t I I I 1 I I t I I I l I I i I I i I I I I I ! I I I I t I I I I l 

n REAL::OPERATOR+() ADD A REAL TO ANOTHER REAL 

///>/////////////////////////////////////////////////////////////// 
// Brief description: add a real to another real 
// Expected: The current real class 
// Result: real reference 

/////////////////////////////////////////////////////////////////// 

real real :: operator + (real i) 

{ /* real + real*/ 

real r = *this; 

return r += i ; } 

/////////////////////////////////////////////////////////////////// 

i REAL: : OPERATOR<() COMPARE TWO REALS FOR LESS THAN 

/////////////////////////////////////////////////////////////////// 
// Brief description: compare two reals for less than 
// Expected: The current real class 
// Result: real reference 

/////////////////////////////////////////////////////////////////// 
int real :: operator < (real i) 
{ /* real < real*/ 

return lss (this, &i);} 

/////////////////////////////////////////////////////////////////// 
„ REAL::OPERATOR<=() COMPARE TWO REALS FOR LESS THAN OR EQUAL 

/////////////////////////////////////////////////////////////////// 
// Brief description: compare two reals for less than or equal 
// Expected: The current real class 
// Result: real reference 

//////////////////////;//////////////////////////////////////////// 
int real :: operator <= (real i) 
{ /* real <= real*/ 

return leq (this, &i);} 

/////////////////////////////////////////////////////////////////// 

// REAL::OPERATOR<=() COMPARE TWO REALS FOR GREATER THAN 

/////////////////////////////////////////////////////////////////// 
// Brief description: compare two reals for greater than 
// Expected: The current real class 
// Result: real reference 

/////////////////////////////////////////////////////////////////// 

int real :: operator > (real i) 

{ /* real > real*/ 

return lss (&i, this); /* if 'i' is less then we are greater*/} 

/ / / / ,.,./._,/.._/ . / / // // // // // // // // // // // // // // // // // // // // // // // // // // // // // 

// REAL::OPERATOR>=() COMPARE TWO REALS FOR GREATER THAN OR EQUAL 

/////////////////////////////////////////////////////////////////// 
// Brief description: compare two reals for greater than or equal 
// Expected: The current real class 
// Result: real reference 

/////////////////////////////////////////////////////////////////// 

int real :: operator >= (real i) 

{ /* real >= real*/ 

if (lss (this, &i)) /* we are less*/ 

return 0; 

return 1? / * we must be >=*/} 

/////////////////////////////////////////////////////////////////// 

n REAL::OPERATOR%() THE MODULUS REMAINDER OF A REAL AND AN INTEGER 

/////////////////////////////////////////////////////////////////// 
// Brief description: the modulus remainder of a real and an integer 
// Expected: The current real class 
// Result: real reference 

/////////////////////////////////////////////////////////////////// 
real real :: operator % (int i) 
{ real a = (*this / (real) i) ; 

uint64 scale [J = { SCALE, 0 }; 

divl28 (this, (real *) scale, this, 0) ; 

real r = *this - (real) (a * (int64) i) ; 

return r; } 

/////////////////////////////////////////////////////////////////// 
h REAL::OPERA TOR*() MULTIPLY TWO REALS 

/////////////// //////////////////////////////////////////////////// 
// Brief description: multiply two reals 
// Expected: The current real class 
// Result: real reference 

/////////////////////////////////////////////////////////////////// 
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real real :: operator * (real m) 
{ real r; 

mulr (this, &m, &r) ; 

return r; } 

///////////////////////////////////////////////////////// ////// ft i 1 
f RE AL: :OP ERATOR&() A BIT AND OF A REAL AND AN INTEGER 

/////////////////////////////////////////// / / ////////////////////// 
// Brief description: a bit and of a real and an integer 
// Expected: The current real class 
// Result: real reference 

////////////////////////////////////////// / iii / // / / / / / // // // // // // / 
int real :: operator & (int m) 

{ return (int) (ip & m) ; /* bitwise AND on integer portion only*/} 

///////////////////////////////////////////////////////// ' / / / / / / / / / 
« REA L::OPERAT OR/( ) DIVIDE TWO REALS 

//////////////////////////////////////////////////////////////'//// 
// Brief description: divide two reals 
// Expected: The current real class 
// Result: real reference 

//////////////////////////////////////// / / / // // / / / / iii i i i i / it ////// 
real real :: operator / (real m) 
{ real r; 

divr (this, &m, &r) ; 

return r; } 

/////////////////////////// / //////////// ////ii/ i / / / // // i / / i i / i i i i i i 
u RE AL::OPERATOR*=() MULTIPLY TWO REALS AND ASSIGN 

~7 7 7 7 7 7 7 7 7 7 7 7 7**7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 

/ / Brief description: multiply two reals and assign 
// Expected: The current real class 
// Result: real reference 

/////////////////////////////////// / ////////////////////// t /////// / 
real real :: operator *= (real r) 
{ return *this = (*this * r) ? } 

i i t i i i i / i i i / / / / i i i i i i t i i i i i i t i i i i i i / i i / i t / / / i i t i i i i i iii i i i i i i i i i / / t 

u RE AL: :OPERATOR/=() DIVIDE TWO REALS AND ASSIGN 

//////// i i // i / i i i / / i i i / i i i / i i i / i i i ////// / // i // i i / / / i i / i / i i i i i i / i i / i 
II Brief description: divide two reals and assign 
// Expected: The current real class 
// Result: real reference 

/////////////////////////////////////////////////////////////////// 
real real :: operator /= (real r) 
{ return *this = (*this / r) ; } 

/////////////////////////////////////////////////////////////////// 
// CMPR() COMPARE TWO REALS FOR EQUALITY 

////////////////////////////// / / // // // / ////// ////// / // // // // / / / / / / / 
II Expected: The address of two reals 
// Result: -1 if a < b 
// 0 if a == b 

// 1 if a > b 

/////////////////////////////////////////////////////////////////// 
int cmpr (real * a, real * b) 
{ if (*a == *b) 
return 0; 
if (*a < *b> 

return -1? 
return 1 ; } 

/ //////// i i / / i / t i i i i i i i i / i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i t i i i / i ii t i 



n REAL::OPERATOR==() COMPARE TWO REALS FOR EQUALITY 



i i / / i // i // / / i i / i i / i / i i i i // i i i i i i i i i i i / i / i / i i 
If Expected: The current real class 
// Result: real reference 

i i i / i / i i / i / i i i i i / / i i i i / i / i / i i i i / i i i / i i i / t i / i 

int real :: operator == (real r) { return ip == r.ip && fp == r.fp;} 

i / i i i i i i / / / / / / // / / i / / i i i i i i / i i i / i i i / i / i i i i i i 



////// 



i i i i i i 



ii/ii/ 



iiiii/ii/ 



i i i i i i / i i 
i / / i i i i / i 



i i i 
/ / i 



n REAL::OPERATOR!=() COMPARE TWO REALS FOR INEQUALITY 



i i i i i i i i i 



i i i / i i i i i 



i / i i i i / i i 



i / / i / i / i / i i / i i / i i / / i / i / / i i i / i i i i i i i i i i i i i i i i 
II Expected: The current real class 
// Result: real reference 

//////////////////////////////////////////// 
int real :: operator != (real r) 
{ return ! operator == (r) ; } 

I I I I I i i / / i // i i / i / / t / / / // t i // i i / I i I I I I I i i i / i i 



i i i i i i 



i i i i i / 



i i / i i i 



REAL::OPERATOR!() LOGICAL NOT OF A REAL (IS ZERO) 



/ i i i i i i i / i / i i // i / / / / i / / i / i i / i i i i i i i i i /// / / / i 
II Expected: The current real class 
// Result: real reference 

i i i / i i / / i / i i i / i i i / i i i i i i i / i t i / / / / / / / i i / / / / / i 
int real :: operator ! () 
{ return I ip && I f p ; } 

//////////////////////////////////////////// 



i i i i i i 



i i i i i i 



i i i i / i 



i i i i i i i i i 



i i i i i i i i i i i i / 



i i i i i i i i i i i / i 



Appendix, Page 26 of 43 



QIPLG Attorney Docket No:500.0028-30us 



m REAL: :OPERATOR»=() RIGHT BIT SHIF T OF A REAL AND ASSIGN 

//////////////////////// ///////////'ft'/'''''' 11 ' 1 f 1 f / ' 1 / f 1 / 1 1 ' ' 1 1 1 
II This function has not been implemented (it is not used) 
// Expected: The current real class 
// Result: real reference 

I , f ( t t t f ! I I 1 f l ! t i I I i I I I t I I I 1 I t t I I I I t t I l I ' I I < ' I t I < 1 1 1 1 1 1 1 1 1 ' 1 ' 1 ' f 1 f ' 

real real :: operator >>= (int) 
{ return *this; } 

/////////////////////// 



/////////////////// / i iii iii iii iii i i i i i i i i i i i 



a REAL::OPERATOR»=() LEFT BIT SHIFT OF A REAL AND ASSIGN 

i , , , , i i , , i , , i i , , i i i i i i i i i i i i i i i i i i i i i i i i i i i i i { i i f ' i i i i ' < i i ' ' 1 1 1 > 1 ' ' 
II This function has not been implemented (it is not used) 
// Expected: The current real class 
// Result: real reference 

y t i i t / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / * / i t 1 1 1 f 1 1 1 1 ' / ' / 
real real :: operator <<= (int) 

{ return^ this,^ t f t t f f t , , t , f f , , ( f , , , t t t , , , , t t i i i t t t i t i i t i i i i » i t i i i i i t t i 



n REAL:: OP ERATO R»=() LEFT BIT SHI FT OF A REAL 

( 1 ! I { i i l ! ! t I I I I I I I ! f i I i I J I I I I I I 1 i I I i I I I I I I I I I I I I ' ' I I t I I I I I I I 1 I I l I l I 

II This function has not been implemented (it is not used) 
// Expected: The current real class 

// Result: real reference ,,,,,,,,,,, 
, f f t t i t t i i t i i t i i i i i t i i i t t i i t i t i i i t i i t t i t t t i i i ( i t i t t t t i t t t i i * t i f t i t t 

real real :: operator << (int i) 

{ real r = *this << i; 

return r; } 

i i i i i i t i i i i i t i i i i i i i i i i 



i i i i i i i i i i i i i i i i i i i i i i i i f t i i iii iii i i i i i i i iii 



n REAL::OPERATOR»={) RI GHT BIT SHIFT OF A REAL 

t I f I I I ! f I I t i I I l I l I I I I I l I I l I I I I I I i I I I I i I I ! t I I I I 1 I I I I I I I 1 I I I I I I i ! I l i I 

3 „ // This function has not been implemented (it is not used) 

til // Expected: The current real class 

O // Result: real reference 

ZZ / / t / / / i / / f / / / / t i f / t i / / t f / t i / / / t / i i f / / f / / t f t t ( * t t t * ' * ( t t t I / t / / / I f f f i 

Q real real :: operator >> (int i) 

if! { real r « *this >> i; 
^ return r;} 

%kS jl for rtoaO ascii digit table 

yQ #define TABLE_S C ALE 10000 

:S #define TABLE_DIGITS 4 /* number of »0's*/ 

^ static char ntab [TABLE SCALE] [TABLE DIGITS] j 

i i i i i i i I i i i i i i i / i i i i i i i ~i t i i i i i i i i i i i i i i i i i i iii ii i i i t i i i i i i i i i i i i i i I 



// MATHJNIT 

t f I I I ~l I I I I I ! I t I l I I l I I I I i I I l I t l l I ' I l I I * I * * I * I I I * i i ' 1 f > 1 1 > > 1 1 ' ' ' 1 1 1 1 

f% If Brief description: Initialize numeric table with ascii values 
rf. II The table is initialized with the least significant bytes 
W // first. For example, entry ntab[l] would be 

yk II ntab[l] [0] = '1'; // least significant digit first 

// ntab[l] [1] = '0'; 

^ // ntab [13 [2] - '0' ; 

Q // ntab[l] [3] = '0«; 

fll // and entry number ntab [1678] would be 

55# // ntab [1678] [0] = '1'; // least significant digit first 

// ntab [1678] [1] = '6' ; 

// ntab [1678] [2] = »7' ; 

// ntab [1678] [3] = '8' ; 

// Result: void 

////////////////////////////////////////////////////////////////''' 
void routine math_init () 
{ char ascii_buf [32] ; 

char format_string [32] ; 

int i, n; 

sprintf (format_string, "%%0%dd", TABLE_DIGITS) ; /* build format string*/ 
for (i = 0; i < TABLE_SCALE; i++) 

{ sprintf (ascii_buf, format_string, i) ; /* get ascii digits*/ 

for (n = 0; n < TABLE_DIGITS ; n++) /* move those digits into table*/ 

ntab[i] [n] = ascii_buf [TABLE_DIGITS - (n + 1)]; /* least significant first*/ 

J - • i i i i i i i i i i i i i i i i i i i i i t i i t i i i i i i i i i i i i i i i i i t f i f f i i ' i i i i i ' f i i i 



i i i i i i i 



» GET_ASCII_FRACTIONAL_DIGITS 

//////// f i / i i i / i i / i i i i i i i i i i i i i iii i iii i iii i i i ii i i i i i / i i i i i / i i i i i i i i 
II Brief description: Get ascii fractional digits. 
// Expected: 

// a = address a character array for ascii digits 

// i = index into above character array 

// n = binary integer to convert 

// p = scale fractional digits 
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// digits = number of precision digits 

// Result: i new index into a 

//////////////////////////////////////////////////////////////////' 
inline int routine get_ascii_f ractional_digits (char *a, int i, int n, 

int f ract_digits, int digits) 

{ int table_digit, index ; 

if (n < 0) /* calculate -1*/ 

n - -n; /* as 1*/ 

while (f ract_digits > 0) 

{ if (n < TABLE_S CALE ) 

{ /* already less no need for a divide*/ 

index = (int) n; /* use the number as the index*/ 
n = 0; /* no more number*/ } 

else 

{ index = n % TABLE_SCALE ; /* remainder*/ 

n /= TABLE_SCALE; /* new n*/ } 
for (table_digit = 0; table_digit < TABLEJDIGITS ; 
table_digit++, --f ract_digits) 
if (fract_digits > 0 && fract_digits <= digits) 
a[i++] = ntab [index] [table_digit] ; } 

return i;} /* get ascii fractional digits*/ 

/////////////////////////////////////////////////////////////////^/ 



/; GET ASCII WHOLE DIGITS 



/////////////////////////////////////////////////////////////////// 
// Brief description: Get ascii whole integer digits. 
// Expected: 

y, // a = address a character array for ascii digits 

L,, // i = index into above character array 

%»f // n = binary 64 -bit integer to convert 

O // Result: i new index into a 

*~Z. /////////////////////////////////////////////////////////////////// 

£y inline int routine get__ascii_whole_digits (char *a, int i, real * _n) 

r.=l ( ^ nt index, table_digit, skip_leading; 

'% int 64 n = _n->ip; /* integer portion*/ 

d do 

jPJI { if (n < TABLE_S CALE ) 

fZ { /* already less no need for a divide*/ 

index = (int) n; /* use the number as the index*/ 
a n = 0; /* no more number*/ } 

d»j else 

7 5r ; { index = (int) (n % TABLE_SCALE) ; /* remainder*/ 

W n /= TABLE_SCALE; /* divide by TABLE_SCALE*/ } 

jy. a[i++] = ntab [index] [0]; /* always at least 1 character*/ 

hi if (n) 

= 3 { /* more to come*/ 

O for (table_digit = 1; table_digit < TABLE_DIGITS ; 

ff| table_digit++) 

a[i++] = ntab [index] [table_digit] ; /* convert all digits*/ } 

else 

{ /* skip leading '0's on last divide*/ 

for (skip_leading = TABLE_DIGITS - 1; skip_leading; 
--skip_leading) 
if (ntab [index] [skip_leading] != '0') 

{ for (table_digit = 1; table_digit <= skip_leading; 

table_digit++) 
a[i++] = ntab [index] [table_digit] ; 
break; 

} 

} } 

while (n) ; 

return i;} /* get ascii whole digits*/ 

/////////////////////////////////////////////////////////////////// 



„ RTOA REAL TO ASCII 

/////////////////////////////////////////////////////////////////// 

// Brief description: Convert a real for ASCII output. 

// This routine is almost identical to rtoa. . .except that it fills 

// in real info with information about the converted number. Information 

// includes stuff like: whole_digits , f ract_digits , neg or not, ... 

// Expected: 

// in = address of 64-bit real 

// _desc = address of output descriptor 

// fdigits = number of fractional digits to display 

// flags = (2 = implied decimal point) ??fix this?? 

// Result: real_info is filled 
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z / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / ' ' ' ' ' ' 1 ' ' ' / ' 1 ' 1 ' ' 1 ' ' 7 1 1 ' ' 1 
void routine rtoa (real * in, void *_desc, int fdigits, realinfo * rinfo, 

int flags) 
{ DESC_S *desc = (DESC_S *) _desc; 

char *r = desc->dsc$a_pointer ; 

int len = desc->dsc$w_length; 

int skip_trailing, i = 0; 

char buf [100] ; 

char *a = buf; 

real n; 

int f; 

int s = 0; 

int p; 

int is_implied_dec = (flags & 2) ; 

fill (sizeof (realinfo), rinfo, 0); /* zero all of the realinfo structure*/ 
n = *in; 

if (is_negative_real (&n) ) 
{ negate_real (&n) ; 

b = 1; } 
rinfo- >neg = s; 

f = n.fp; /* get fractional part*/ 

/* * we build the string from left to right 

* in reverse order fraction. whole digits*/ 
if (fdigits || f != 0) 

{ /* fraction first right to left*/ 

H= for (p = SCALE_DIGITS ; p < fdigits ; p++) 

a[i++] = '0'; /* fill with '0' digits first*/ 

!jj i = get_ascii_fractional_digits (a, i, f, SCALE_DIGITS , fdigits) ; 

%J for (skip_tr ailing = 0; 

f ,f| skip_trailing < i && a [skip_trailing] == '0»; skip_trailing++) ; /* skip ■ O's*/ 

rinfo- >fractdigits = i - skip_trailing ; /* number of fractional digits ignoring trailing '0's*/ 

if ( i is implied dec) 

a[i++] = '.'? /* add a decimal point*/ } 

A if (n.ip != 0) 

"f^ { /* get the whole ascii digits (all of them)*/ 

yj int new_i = get__ascii_whole_digits (a, i, &n) ; 

B rinfo- >wholedigits = new_i - i; /* number of whole digits*/ 

rs» i = new_i; } 

^ if (s) 

yj if (i > 0) 

iU a[i++] = 

L. rinfo->rlen = i; /* length of number (relative to 1}*/ 

SU a[i] = '\0'; // for easy debug 

for (p = 0; p < len - i; p++) // blank fill result 
m r[p] = ' '; 

for (; p < len; p++) 

{ // fill result with answer 

r[p] = a[--i] ? // from reverse order 
}} // rtoa 

char* routine rtoa (real *n, int d) 
{ static char buf [28] ; 
DESC_S desc [1] ; 
realinfo info; 
desc->dsc$a_j)ointer - buf; 
desc->dsc$w_length = sizeof (buf) ; 
rtoa (n, desc, d, fcinfo, 0) ? /* real to ascii*/ 
return buf ; } 

////////////////////////////////////////////////// ft ''I' i if t i I L ' 1 1 1 



„ NAT_CNVOUT CONVERT A REAL FOR A SCII OUTPUT 

///////////////////////// /////// / / / f i f f f f f f f f f t f ' f f ' f f i 1 1 / 1 1 1 1 ' ' ' 1 1 
II Brief description: Convert a real for ASCII output. 
// Expected: 

// val = address of 64 -bit real 

// result = address of ASCII buffer 

// digits = number of digits to display 

// Result: status from math cnvout 

/,iff i i /////// i t / / //">///////////// / / /////// f / / / / / / / / f f f f f f f i f ' 1 f ' 1 ' 
int routine math_cnvout (real * val, void *result, int digits) 
{ DESC_S *desc = (DESC_S *) result; 

char *dest = desc- >dsc$a_j?ointer ; 

int rlen = desc->dsc$w_length; 

char tbuf[48] ; //?? why 48 

int len = 38; //?? why 38?? 
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int offset; 

DESC_S tmp (tbuf, len); 
real info rinfo; 

fill (rlen, dest, 1 '); // blank-fill result 
rtoa (val, &tmp, digits, fcrinfo, 0) ; 

offset = len - rlen; //?? what if RESULT buffer is bigger than ours? 
if (offset < 0) 

offset = 0; 
copy (rlen, &tbuf [of f set] , dest) ; 

, i f e / t y r / n / T / R V E /l ////////////;///////////////////// LA L LU L-L ' ' t-L I LA L L > L L L L L 
1 u NEW CNVOUT CONVERT A REAL FOR ASC II OUTPUT 

1 7-7-7-7 7 ///////////////////////////// 11 'I 1 1 ? 1 1 ' 1 1 7 1 1 ' 1 ' ' 1 1 1 ' ' 1 1 1 ' ' 1 ' 1 1 
II Brief description: Convert a real for ASCII output. 
// Expected: 

// val - address of 64-bit real 

// result = address of ASCII buffer 

// digits = number of digits to display 

// flags = (2 = implied decimal point) ??fix this to be a symbol?? 

// Result: realinfo structure is filled and returned ,,,,,,,,,,,,,,,,/,, 
/////////////////////////////// / / // / 111 1 t / j t / 1 1 1 1 t 1 1 1 1 1 1 1 1 t 1 1 1 1 1 1 * 1 
static int routine new_cnvout (real * val, void *result, int digits, 

realinfo * rinfo, int flags, int maxJLen) 
{ DESC_S *desc = (DESC_S *) result; 

char *dest = desc->dsc$a_pointer ; 

int rlen = desc- >dsc$w_length; 

char tbuf[48]; //?? why 48 

int len = 38; //?? why 38?? 

M= int offset; 

DESC_S tmp (tbuf, len) ; 

fill (rlen, dest, ' '); // blank-fill result 
O rtoa (val, &tmp, digits, rinfo, flags) ; 

if* offset = len - rlen; //?? what if RESULT buffer is bigger than ours? 

^ if (offset < 0) 

5 =?y offset - 0; 

i&g if (max__len < rinfo- >rlen) { fill (rlen, dest, »*'); 

rr% return FALSE; 

lH copy (rlen, &tbuf [of f set] , dest); 

s return TRUE;} // new_cnvout 

f*\ const int cvtnumlen = 32; 

r*"' static int ten_power[] = { 1, I* 0*1 

yy 10, /* 1*/ 

U 100, /* 2*1 

hi 1000, /* 3*1 

10000, /* 4*/ 

O 100000, /* 5*/ 

m 1000000, /* 6*/ 

= " ::;r: 10000000, /* 7*/ 

100000000, /* 8*/ 

1000000000 /* 9*1}; 

//// ////////////// 1 / 1 1 /// s / 1 1 1 1 1 1 1 s 1 1 1 1 1 1 / ' 1 1 1 ' s / 1 ' / 1 ' ' ' 1 1 1 ' i i ' J i i 1 



SCALESREAL SCALE A REA L TO A NUMBER OF FRACTIONAL DIGITS 

/,/////////////////// 1 1 1 1 / 1 1 1 s 1 1 1 1 / 1 1 1 1 1 1 1 ' 1 1 1 ' 1 ' 1 ' ' 1 ' ' 1 1 1 i ' 1 1 ' i 1 1 ' 
II Brief description: Scale a real to a number of fractional digits 
// Expected: 

// out = address of 64 -bit real 

// in = address of 64-bit real 

// frac = number of fractional digits (scale to) 

// Result: void 

// Example: 

// scale 87.36 by 1 = 87.4 
// scale 87.36 by 0 = 87 
// scale 87.36 by -1 = 90 
// scale 87.36 by -2 = 100 

/ , / // // // // // // // // // // // // // / ////// / / / / / / / / / / / ft 1 ft t ft t ti f t 1 / / 11 / / 
void old$scale$real (real * out, real * in, int frac) 
{ int64 scale [] = { SCALE, 0 }; 

int64 one[] = { 1, 0 }; 

int neg = 0; 

if (is_negative_real (in)) 
{ negate_real (in) ; 

neg = 1; } 
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if (frac < -SCALE_DIGITS) 

frac = -SCALE_DIGITS; 
if (frac > SCALE_DIGITS) // our maximum precision 

frac = SCALEJDIGITS; 
/* * here we calculate the scale*/ 
if (frac < 0) 

{ /* -i == o.l*/ 

int64 pow[] = { ten_power [-frac] , 0 }; 
mull28 (scale, pow, scale) ; } 

else 

{ int64 pow[] = { ten_j>ower [f rac] , 0 }; 

divl28 (scale, pow, scale, 0) ; } 
if ( ! is_negative (scale, 2)) 

{ /* not negative*/ 

int64 n[] = { 0, 0 }; /* whole 'n**/ 
int64 r[] = { 0, 0 }; /* remainder 'r'*/ 
int64 a[] = { 0, 0 }; 
int64 b[] = { 0, 0 }; 
a[0] = in->ip; 
scaleup (a) ; 
b[0] = in->fp; 

add2 (a, b) ; /* join ip & fraction*/ 

divl28 (a, scale, n, r) ; 
shiftdown (scale, 2) ,- 
if (leq (scale, r, 2) ) 

add2 (n, one) ; 
shiftup (scale, 2) ; 
mull28 (n, scale, n) ; 
scale [0] = SCALE; 
scale [1] = 0; 

divl28 (n, scale, a, b) ; /* split the 128-bit integer*/ 
out->ip = a[0]; /* to integral part*/ 
out->fp = tint) b[0]; /* fractional part*/ } 
if (neg) 

/ / / ? 6 / 9 f ^T 1 ?^ 1 / \°^) '\ I / I I J l I I l / t I ! I f i i I t i I I I i i t i ' i ' i i J i * ' i t i > i 1 ' ' 1 ! ' ' 1 ' 



ROUND$REAL ROUND A REAL TO A NUMB ER OF FRACTIONAL DIGITS 

/ 7 / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / 7 / t i / t / f i / i t / / / / / / / / 7 / / / / 



Q // Brief description: Round a real to a number of fractional digits 
r v = // Expected: 
w // out = address of 64-bit real 

hh If in = address of 64 -bit real 

m II frac = number of fractional digits (scale to) 

]JZ II Result: void 
W // Example: 

ffj // Round 176.357143 to 2 fractional digits = 176.36 
// round 176.357143 to -1 fractional digits = 180 
// scale 87.36 by 1 = 87.4 
// scale 87.36 by 0 = 87 
// scale 87.36 by -1 = 90 
// scale 87.36 by -2 = 100 

w i , i ( i i { t i t i t / i t i t i t i t t i i i i i i i t i i i i i t i i t t i i i f i i ' 1 1 1 11 'I ' 1 ' 1 1 1 ' ' ' 1 1 
void routine round$real (real *out, real *_in, int frac) 
{ real rin = *_in; 

real *in = &rin; 

int scale = SCALE ; 

int rem; /* remainder*/ 

int neg = 0; 

if (is_negative_real (in)) 
{ negate_real (in) ; 

neg =1? } 
if (frac < -SCALE_DIGITS) 

frac = -SCALE_DIGITS; 
if (frac > SCALE_DIGITS) /* our maximum precision*/ 

frac = SCALE_DIGITS; 
/* * positive 'frac' then only scale the »fp' part*/ 
if (frac == 0) 

{ out->ip = in->ip; /* copy over ip part*/ 

if ((SCALE >> 1) <= in->fp) /* ■ f p ' part greater than a half*/ 

out->ip += 1; /* carry one*/ 

out->fp =0; /* no decimal digits*/ } 
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else if (frac > 0) 

{ out->ip = in->ip; /* copy over ip part*/ 

scale /= ten_j>ower [frac] ; 
rem = in->fp % scale; /* get remainder*/ 
out->fp = in->fp / scale; /* scale down*/ 
if {(scale >> 1) <= rem) /* is there a carry?*/ 

out->fp += 1; /* carry one*/ 

out->fp *= scale; /* and up again (to get rid of lower digits)*/ 

if {out->fp >= SCALE) 
{ /* overflow ?*/ 

OUt->fp -= SCALE; 
out->ip++; /* carry one*/ 

} } 

else 

{ /* scale the 'ip' part*/ 

int64 ip = in->ip; 

if ( (SCALE >> 1) <= in->fp) /* ' f p 1 part greater than a half*/ 

ip += l; /* carry one*/ 

out->fp ==0; /* zero result fp part*/ 

scale = ten_jpower [-frac] ; /* a ten_jpower*/ 
rem = (int) (ip % scale); /* get remainder*/ 
out->ip = ip / scale; /* scale down*/ 
if ( (scale >> 1) <= rem) /* is there a carry?*/ 

out->ip 1; /* carry one*/ 

out->ip *= scale; /* and up again (to get rid of lower digits)*/ } 

if (neg) 

■ / / ^/^/V?* 1 ; /^V / \ ;//////////;///////////////'' LJ. ft LJ LJ > ' ' ' I L ' ' 1 1 ' 1 ' 



«| 0 SCALE$REAL SCALE A REAL TO A NUMBER OF FRACTIONAL DIGITS 



//////////////////////////////////////////// / ^ / ^ > ~f f ' 1 I' ! ' ' < ' ' / ' 1 1 ' 1 
O // Brief description: Scale a real to a number of fractional digits 
,fx (f Expected: 

% 5 // out = address of 64-bit real 

ypj // in = address of 64 -bit real 

,r= // frac = number of fractional digits (scale to) 

O / / R ? S / U /V / ////////////////////////// * I'l i ii 1 > i iii ' 1 1 > 1 1 7 ' 7 ' ' 1 1 1 1 ' 

void scale$real (real * out, real * in, int frac) 
y? { round$real ((real *) out, (real *) in, -frac);} 
s char *dotptr; 
f 8 ! char *cvtnumptr; 

static int zeros, leader, frac, sci, period, neg, pad, leadzero, e; ,,,,,,,,,,,,,,,,,,,, 
= i /////// /////////////////////////// / / / / / / / / LJ. ii i 1 i i i i i i i i ' i i i i i i i ' 1 1 



n FL TO STRING FRAC . 

////////////////////////////////////////// //////// ii i ii i i i i 1 i i i ' i i i 

If Brief description: Convert a float to a string 

// Expected: val - value to be converted to a string 

// outptr - an address of conversion. 

// len - maximum length of conversion 

// digits - the number of fractional digits 

// flaqs: 

//////////////////////////////////// / / / / / / / / / / / / ' ' / ' ////// / ' tiii t * t 
int fl_string_output (int *len, char *outptr) ; 

int routine new_f l_string_output (int *len, char *outptr, int flags) ; 

int routine f l_to_string_f rac (int *len, char *outptr, real * val, int f digits, 

int flags) 
{ char tmpbuf [cvtnumlen] ; 

DESC_S tmpdscEl]; 

real rvalue; 

int status; 

realinfo rinfo; 

int is_implied_dec = (flags & 2) ; 
int *s = (int*) &val->ip; 
*outptr = '\0 ' ; 
if (len[0] <= 1) 
return 0; 

set$desc (tmpdsc, cvtnumlen, tmpbuf) ; 
pad = (flags & cvt_pad) != 0; 
rvalue - *val ; 
s = (int*) &val->ip; 

new_cnvout (&rvalue, tmpdsc, fdigits, &rinfo, flags, *len) ; 
if (tmpbuf [0] == '*') 
return 0; 

zeros = fdigits - rinfo. fractdigits ; 
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dotptr = &tmpbuf [cvtnumlen - rinfo.rlen] ; // point us to the first data 
neg = r info. neg; 
if (neg) 

dotptr = &dotptr[l]; // skip sign (??dme?? fix this later) 
sci = FALSE; 

leader = rinfo.wholedigits; 
period = leadzero = 0; 

frac = f digits; // rinf o. f ractdigits ? 

I 111 fix & 2 to be "IMPLIED_DECIMAL" 
if (is_implied_dec) 

status = new_f l_string_output (len, outptr, flags) ; 
else 

status = f l_string_output (len, outptr) ; 
if ('status) 

fill (len[0] , outptr, »*' ) ; 
return 1;) // fl to string frac 

i ! i i i i t I { t l f i t t I I i t I i i i t i t t~ i~ t i f i t i i > t 1 1 1 > 1 1 ' 1 1 1 1 1 * 1 1 1 ' ' ' ' ' 1 ' 1 ' 1 1 1 

| o CVT OUT D F 

' / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / ' i * t * 1 1 1 1 ' 1 1 1 1 1 1 ' 1 ' 1 1 1 ' 1 
II Brief description: Convert real to floating ASCII output. 
// Expected: 

// fval = address of 64 -bit real 

// dlen = fill length 

// desc = pointer to output descriptor 

// Result: status from fl to string frac ,,,,,,,,,,,,,,//,/ 

int routine cvt_out_d_f (real * fval, int dlen, void *desc) 

{ int status = SUCCESS ; 
N= int len; 

f*| char *text; 

char numbuf [100] ; 
U split$desc (desc, &len # fctext) ; 

lQ round$real (fval, fval, 0) ; 

~~ status = f l_to_string_f rac (&len, numbuf, fval, 0, 0) ; 

if (! status) 
y!j return MSG_CVTERROR ; 

// len -= 1; 
"rZ fill (dlen, text, ' ') ; 

Ui copy (len, numbuf, text + dlen - len) ; 

return status;} // cvt__out_d_f 

/////////////////////////////////// /// L L i i ft i L i 1 !—L-L L i L L L 1 1 ' ' ' 1 / 7 7 ' ' 



« REAL CVT OUT D F 

/;////////////////////// ///////////////////'/////' / / / / / / / / / / / / 



// Brief description: Convert real to floating ASCII output. 
// Expected: 

// fval = address of 64 -bit real 

// dlen = fill length 

// desc = pointer to output descriptor 

// Result: status from fl to string frac 

int routine real_cvt_out_d_f (real * fval, int dlen, void *desc) 
{ int status = SUCCESS; 

int len; 

char *text; 

char numbuf [10 0] ; 

char *p = numbuf; 

split$desc (desc, &len, &text) ; 

round$real (fval, fval, 0) ; 

status = fl_to_string_frac (&len, numbuf, fval, 0, 0) ; 
if (! status) 

return MSG_CVTERROR ; 
fill (dlen, text, ' « ) ; 
if (dlen < len) 

{ p += len - dlen; 

len = dlen; } 
copy (len, p, text + dlen - len) ; 

return status ; } ........ 

{ I I / l I I I I I I I t t I t f I 1 I l I I I I I i t I I f t I t l I t i t I I ! I t f > I I t I I I I t > t I t t I I t i / ) t t 



n NAT CVT_STR_FLT CONVERT ASCII FLOATING S TRING TO REAL 

T7 / fn /////// f / / / t // t /// / // / / / t / t /// t ' ' i //////>// it ii * i i i t i / ' ' * i f > 1 1 
II Brief description: Convert ascii floating string to real 
// Expected: 

// result = address of result 64 -bit real 

// dlen = length of ascii text 
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If text = ascii input text 

, W^/V / ?MM / / ////// / / / / / / / / / / / / / / / ///////////// > " 

int routine math_cvt_str__f It (real * result, int len, char *text) 
{ real n = 0; 
int neg = 0; 
int e = 0; 
int i; 

for (i = 0; i < len text [i] && text [i] <= ' * ; i++) ; /* skip leading*/ 
if (text[i] == '-') 
{ neg = 1; 

} 

for (; i < len && text [i] && text [i] <= ' '; ; /* skip leading*/ 

if (i >= len) 

{ *result = n; 

return 1; } 

if ("(tolower (text[i]) >= '0' && tolower (text [i] ) <= '9')) 
{ // not numeric 

switch (tolower (text[i])) 
{ case 1 e • : 

case 'x' : 
case ' + ' : 
case ' . • : 
break; 
default : 

return 0; 

} } 

D /* * get the integral part first*/ 

S for (; i < len && text [i] >= '0' && text[i] <= '9'; i++) 

W { n.ip *= 10; 

; yQ n.ip += textti] & OxOf; } 

if (i < len && tolower (text[i3) == 'e') 
{ e = 0; 

for (i++; i < len && text Ei] >= '0' && text [i 3 <= '9*; 
{ e *= 10; 

e += text [i] & OxOf ; } 
while (e--) 

n.ip *= 10; } 
if (i < len && textti] == ' 
{ int f » 0; 

W int p = ++i; 

jU if (i < len ScSc (textti] < '0' || text[i] > '9')) 

r: return 0; 

*if for (; 

O i < len && i < p + SCALE_DIGITS && textti] >= ' 0' 

fl§5 && textti] <= '9'; i++) 

si# { f *o 10; 

f += text [i] & OxOf; } 
/* do not round ascii input*/ 

/* if (i == p+6 && textti] >= '5' && textti] <= '9')*/ 
/* f += 1; */ 
for (; i < p + SCALE_DIGITS; i++) 

f *= 10; 
n.fp = f; } 
if (i < len && tolower (textti]) == 'e') 
{ e = 0; 

for (i++; i < len && textti] >= '0' textti] <= '9'; i++) 
{ e *= 10; 

e +« textti] & OxOf; } 
while (e--) 
{ n.fp *= 10; 

n.ip *= 10; 
if (n.fp >= SCALE) 

{ n.fp -= SCALE; 

n. ip++; 

} 

} } 

if (is_negative_real (&n) ) 

{ runtime$error (MSG_FLOATOVF) ; /* overflow*/ } 

if (neg) 

negate_real (&n) ; 
*result = n; 
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return 1 ? } 
int routine ator(real *a, char *b) 

I j NAT_CyT_STR_FLT IMPDEC CONVERT ASCII FLOATING STRING TO REAL 

// Brief description: Convert ascii floating string to real 
// 

// Expected: 

// result = address of result 64-bit real 

// dlen - length of ascii text 

// text = ascii input text 

, , f' a w s / / / //////// / / /////////////// / / / ' / ' .' ' ''/''' " " " 11 ' ' 1 ' 

int routine math_cvt_str_f lt_impdec ( real * result, int len, char *text, mt f digits) 
{ real n = 0; 
int neg = 0; 

int i ; „ . , , . ^ , 

for (i = 0; i < len && text[i] && text[i] <= ' '; i++) ; /* skip leading*/ 

if (text[ij == 1 -■) 
{ neg = 1; 

i++; } 

for (; i < len text [i] && text[i] <= ' * ; i++) ; /* skip leading*/ 
if (i >= len) 

{ *result = n; 

return 1; } 
/* * get the integral part first*/ 

for (; i < len-f digits && text[i] >= '0 r && text[i] <= '9'; 
Lfe { n.ip *= 10; 

n.ip += text[i] & OxOf; } 
S int f = 0; 

O int p = i; 

if* if (i < len && (text[i] < '0' || text [i] > '9')) 

'% return 0; 

^ for (; i < len && i < p + SCALE_DIGITS && text [i] >= '0' 

hf\ && text[i] <= '9'; i++) 

S { f *= 10 ? 

"rt f += text[i] & OxOf; } 

HI for {; i < p + SCALE_DIGITS; i++) 

s f * = 10 '* 

n.fp = f; 
W if (is_negative_real (&n) ) 

hj { runtime$error { MSG_FLOATOVF ) ; /* overflow*/ } 

u if {neg) 

L,. negate_real (&n) ; 

IU *result = n; 

O , , f^/V/ } , / / /////////////////// /////// UL ' ' UL > ' ' ' < < 1 * L 1 '- / I ' 1 1 f 1 { 1 1 1 1 1 

% U DTOR CO NVERT A DOUBLE TO REAL ,,,,,,,,,,,,,,, t i i 

1 //////// 7 ////////////////////////////'' ft ' 1 i ! ' ' 1 1 ' ' ' / 1 1 1 / ' ' 1 1 1 ' f 1 1 ' 

II Brief description: Convert a double to a 12 8 -bit real 

// Expected: 

// a address of input double. 

// b address of output 128-bit real 

void 

routine dtor (double *a, real * b) 
{ double d = *a; 

double my_d; 

int64 ip; 

int fp; 

if (d > LARGEST_NUMBER | | d < - LARGE ST_NUMBER) 

runtime$error (MSG_FLOATOVF) ; 
// handle rounding now by adding in 5*10 A -8 
if (d < 0) 

my_d = d - 0.000000005; 
else 

my_d = d + 0.000000005; 

ip = <int64) my_d; 
my_d -= ip; 

fp = (int) (my_d * SCALE); // this is the same code I last sent you 
fp = (fp/10)*10; // truncates at SCALE-1 digits 
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if (ip && fp < 0) /* if ip then*/ 

fp = -fp; /* fp always positive*/ 

b->ip = ip; 
b->fp = fp; 

, / / \ //////////////////////// / ///'''/ L ' ' 11 11 11 1 1 1 1 1 ' 1 ' 1 1 1 1 { ! 1 1 1 ' ' ' ' ' 
u MULR MUL TIPLY TWO REALS ,,,,,,,,///////;///;//;;; 

//////////////////////////// //// rn / / / / / / / / / / / // ' > * < < 1 ' 1 1 ' 1 1 1 1 1 1 111 

II Brief description: Multiply two 128-bit integers. This is used for intermediate 
// arithmetic. 
// Expected: 

// _a = address of real (_a * _b) 

// _b = address of real (_a * Jo) 

// r address of real (result) 

, fVf^l / T^f, /////// / / ///////////// / / / / ////////// " / ' ' ' ' " ' " ' " " ' 

void mulr (real * _a, real * __b, real * _c) 

, { , ggS^ /~^?'/ T~) "'i^ / t i i i i i i i t 1 i i i i i i i i I i i t i i > 1 i { 1 1 1 1 1 1 1 1 ! 1 1 { 1 1 1 1 1 1 1 1 



// CVTSQUADSFLT CONVERT A 64-BIT Q UAPWORD TO REAL ,,,,,,,,,, 

1 7 r-7 / / ////////////////////////////// 11 1 < 1 ' ' 1 > 1 f 1 1 ' 1 ' ' 1 7 1 1 1 1 ' 1 1 ' ' ' ' ' 



////////// / // // // // // /'/' fill f f t t , t < > > 

II Brief description: Convert a 64-bit quadword to a real. 
// Expected: 

// result = address of real 

// p__quad = address of 64 -bit quadword 

ft / JVff; / / st ; a w s / 

y& int routine cvt$quad$flt (real * result, int *p_quad) 
{ result->ip = (*(int64 *) p_quad) ; 

result->fp =0; // Clear fractional portion ++dme++ 



^ / / J^fffVl ////////////////// ///////////// ' 1 11 L ' 1 1 ' t ' ' 1 ! ' ' f ' ' 1 ' 1 ' 1 ' 

1 \ FL TO STRING 



/////////////////////////////////// ' ' rn i j ii / / 7 / / / / / / / / / / / / / / / / / ' ' ' ' 

II Brief description: Convert a float to a string 

// Expected: val - value to be converted to a string 

fl // outptr - an address of conversion. 

TZ II len - maximum length of conversion 

// digits - the number of significant digits 

L i 1 , , / / / / / / f / lag / s / //////////////////////////////'/ { * * < > * 1 < * * * f 1 ' ' ' ' 1 1 1 ( 1 

U int routine fl_to_string (int *len, char *outptr, real * val, int digits, 
Id int flags) 

{ char tmpbuf [cvtnumlen] ; 
^ DESC_S tmpdsctl]; 

|Tj int status; 

.aw ; realinfo rinfo; 

Jrf pad = (flags & cvtjpad) 1=0; 

lU set$desc (tmpdsc, cvtnumlen, tmpbuf) ; 

new_cnvout (val, tmpdsc, digits, &rinfo, 0, *len) ; 

if (tmpbuf [0] == '*') 
return 0; 

zeros = digits - rinfo . fractdigits ; 

dotptr = & tmpbuf [cvtnumlen - rinfo. rlen] ; // point us to the first data 
neg = rinfo.neg; 
if (neg) 

dotptr = &dotptr[l] ; // skip sign (??dme?? fix this later) 
SCi = FALSE; 

leader = rinfo.wholedigits; 
period = leadzero = 0 ; 
frac = rinfo. fractdigits ; 
cvtnumptr = tmpbuf; 

status = f l_string_output (len, outptr) ; 

return status;) // fl to string 

/ / 7 / 7 > i i i \ i / ///////////////////-/-///// //////// LU i < * ' < i > ' L > L ' ' I i L L 1 1 

| \ FL STRING OUTP UT , , , , , 

1 , t t , , , , , , i i t i t i i f i t i i i i i i i t i t i f i > i i > i i i * ' ' i * < 11 1 i 1 1 ' 1 ' ' 1 1 1 / 1 1 ' 1 ' 1 1 1 



II Brief description: Convert a float to a string 
// The result looks like: [- 3 nnnnnnn. f f f f f f f 
// Where [- ] is either a space or '-' 
// nnnn = whole number, f f f f =f raction 
// (notice how the decimal point is included) 
// Expected: val - value to be converted to a string 
// outptr - an address of conversion. 

// len - maximum length of conversion 
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// digits - the number of significant digits 

/ / / / / / / f /^ g / S / 
int routine f l_string_output (int *len, char *outptr) 
{#define $save_lit (ch) \ 
{ \ 

if ( (outptr - outbegin) > maxlen) \ 

return 0; \ 
outptr [0] = ch; \ 
outptr++; \ 

} 

#undef $save 

#define $save(tlen, txt) \ 
{ \ 

if ( {outptr - outbegin) > maxlen) \ 

return 0; \ 
copy(tlen, txt, outptr); \ 
outptr += tlen; \ 

} 

char * outbegin; 
int maxlen; 

memset (outptr, 0, *len) ? /* null fill result*/ 

maxlen = len [0] ; 
outbegin = outptr; 
if (neg) 

$save_lit ('-«) 

else 
if (pad) 

$save_JLit (' '); 
if (leadzero) 

$save_lit ( '0' ) ; 
$save (leader, &dotptr[0]); 
if (frac == 0 && leader == 0 && ! leadzero) 

$save_lit ( ' 0 1 ) ; 
if (frac > 0) 

{ $save_lit ('.'); 

$save (frac, &dotptr[l + leader]); } 
if (pad) 

$save_lit ( ■ ' ) ; 
*len = outptr - outbegin; 
if (maxlen < *len) 

return 0; 

retUrn ,Vl ,,,,,,,,,, 'j ^figSfflL •>>•> ,,,,,,,,,;////,/;/// 



////// 

7 fY* T 7^ / Q / U / T / P / U / T / ///////////////////////////////// ri 7-7 / / / / / / / / n / 7 / 
// Brief description: Convert a float to a string 
// The result looks like: [- 3 nnnnnnnf f f f f f f 
// Where [- 3 is either a space or '-' 
// nnnn = whole number, f f f f =f raction 

// (notice how the decimal point is implied but not included) 

// Expected: val - value to be converted to a string 

// outptr - an address of conversion. 

// len - maximum length of conversion 

// digits - the number of significant digits 

// flags: 2 = implied decimal digits t ,,,,,,, , 

), , i , , t , , , ,, , , , i i i i i i i i i i i t ii t i t i t i i i t i i i i i i i i i t 1 i ( 1 1 1 1 > 1 1 1 1 * 1 1 1 { 1 1 
int routine new__f l_string_output (int *len, char *outptr, int flags) 
{ char * outbegin; 
int maxlen; 

int is_implied_dec = FALSE; 
if (flags & 2) 

is__implied_dec = TRUE; 
maxlen = len [0] ; 
outbegin = outptr; 
if (neg) 

$save_lit ('-') 

else 
if (pad) 

$save_lit ( 1 ' ) ; 
if (leadzero) 

$save_lit CO'); 
$save (leader, fcdotptr [0] ) ; 
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if (frac == 0 && leader == 0 && lleadzero) 

$save_lit CO'); 
if (frac > 0) 

{ if (is_implied_dec) //??dme?? as WHy no ";" on line below?? 

$save (frac, &dotptr [leader] ) // just the decimal digits. . .no . 
else 

$save (frac + 1, &dotptr [leader] ) ; // fract+1 to include the " . " } 

if (pad) 

$save_lit (' '); 
*len = outptr - outbegin; 
if (maxlen < *len + is_implied_dec) 

return 0; 

return 1;} // new_f l_string_output 

/////////////////////////////////////////////////////////////////// 

// FL TO STRING 

/////////////////////////////////////////////////////////////////// 

// Brief description: Convert a float to a string 

// Expected: val - value to be converted to a string 

// outptr - an address of conversion. 

// len - maximum length of conversion 

// digits - the number of significant digits 

// flags: 

/////////////////////////////////////////////////////////////////// 
int routine real_f l_to__string (int *len, char * outptr, real * val, int digits, 

int flags) 
{ char tmpbuf [cvtnumlen] ; 
DESC_S tmpdsc [1] ; 
s a int d = 0/ 

r := int i ; 

Q pad = (flags & cvt_pad) != 0; 
ssat set$desc (tmpdsc, cvtnumlen, tmpbuf) ; 
^ math_cnvout (val, tmpdsc, digits); 
\U if (pad) 

if\ { *outptr++ = ' 1 ; 

■« d = 1; } 

for (i = 0; tmpbuf [i] ; i++, d++) 
Q *outptr++ = tmpbuf [i] ; 

;?g if (pad) 

{ d + +; 
s *outptr++ = ' 1 ; } 

□ *len = d; 
s/jj e = 0; 

frac =0; /* leader - zeros;*/ 

f& sci = 0; 

leader = 0; 

period = leadzero = 0; 
Irj return 1;} // fl_to_string 

jk; : /////////////////////////////////////////////////////////////////// 

*W Q FL TO STRING FRAC 

/////////////////////////////////////////////////////////////////// 

// Brief description: Convert a float to a string 

// Expected: val - value to be converted to a string 

// outptr - an address of conversion. 

// len - maximum length of conversion 

// digits - the number of fractional digits 

// flags: 

// if ( ! f l_to_string_frac (&numlen, buffer, (int*) ^number, 
// f orm->f orm$fraclen, 0) ) 

// return FALSE; 

/////////////////////////////////////////////////////////////////// 
int routine real_f l_to_string_f rac (int *len, char *outptr, real * val, 

int digits, int flags) 

{ char tmpbuf [cvtnumlen] ; 
DESC_S tmpdsc [1] ; 
char *e = outptr + *len; 
char *p = outptr; 

char fill « '0 1 ; /* default*/ 

int i, n; 

if (len[0] <= 1) 

return 0 ; 
pad = (flags & cvt_pad) != 0; 
set$desc (tmpdsc, sizeof (tmpbuf) , tmpbuf) ; 
math_cnvout (val, tmpdsc, digits) ; 



; s s 
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if (pad) 

*p++ = ' 1 / 

for (i = 0; tmpbuf [i] && p < e && tmpbuf [i] != 

*p++ = tmpbuf [i] ; 
if (digits) 

{ if (tmpbuf [i] != ' . ') 

*p++ = 
else 

*p++ = tmpbuf [i++] ; 
for (n = 0; n < digits && tmpbuf [i] && p < e; i++, n++) 

*p++ = tmpbuf [i] ; 
for (; n < digits && p < e; n++) 

*p++ = ' 0 1 ; 
fill = ' ■; } 

else 

{ if (tmpbuf [i] < = « . ') 

{ -e; 

fill = ' ' ; } 
for (; tmpbuf [i] && p < e; i++) 
*p++ = tmpbuf [i] ; } 

if (pad) 

*P = ' 
if (p < e) 

{ /* move to end of string*/ 

while (p > outptr) 

*--e = *--p; 
while (e > outptr) 

*--e = fill; } 

return 1;} // f l_to_string_f rac 

/ / //////////// / // t /////// / t ///////////////////////////// * ///////// / 

h FL STRING OUTPUT 

/////////////////////////////////////////////////////////////////// 

// Brief description.* Convert a float to a string 

// Expected: val - value to be converted to a string 

// outptr - an address of conversion. 

// len maximum length of conversion 

// digits - the number of significant digits 

// flags: 

/////////////////////////////////////////////////////////////////// 
int routine real_f l_string_output (int *len, char *outptr) 
{ char *outbegin; 
int maxlen; 
maxlen = len [0] ; 
outbegin = outptr; 
if (neg) 

$save_lit ( 1 - ' ) 
else 
if (pad) 

$save_lit ( ' 1 ) ; 
if (leadzero) 

$save_lit CO'); 
$save (leader, &dotptr[l]); 
if (frac == 0 && leader == 0 && < leadzero) 
$save_lit CO') 
else 
if {frac > 0) 

{ $save_lit (' . ') ; 

period = 0; 
if (isci) 

while (e < 0) 
{ e++; 
$save_lit CO')} 
$save (frac, &dotptr[l + leader]); } 
else if (isci) 

{ while (e - leader > 1) 

{ e~; 
$save_lit ( ' 0 1 ) } } 
if (period) 

$save_lit ('.'); 
if (sci) 

$save (4, &cvtnumptr [cvtnumlen - 4] ) ; 
if (pad) 
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$save_lit ( ' 1 ) ; 
*len = outptr - outbegin; 

return 1;} // f l_string_output 

# include "tmath.h" 

i i i i i i i i t t I l I I I t ! I I l I I I i I I I I i I I I l I i f ( I f / ( t t / f t t / t t ! I / i t I I I i f I I I I I I I 

it NAT INT SIN 

/////////////////////////////////////////////////// //////////////// 
// Brief description: Returns the SINE of an angle specified in radians 
// Expected: in - passed real an angle 
// Result: out - sine of given angle 

/////////////////////////////////////////////////////////////////// 
double mth$dsin {double *) ; 

c__routine int math_sin (double *in, double *out) 
{ double a; 
real b; 

a = getdouble ((real *) in); 
b = mth$dsin (&a) ; 
* (real *) out - b; 

return 1;} // math_sin 

/////////////////////////////////////////////////////////////////// 

v NAT INT COS 

/////////////////////////////////////////////////////////////////// 
// Brief description: Returns a cosine of an angle. 
// Expected: in - passed real angle 
// Result: out - cosine ofd an angle. 

//////////////////////////////////////////////✓//////////////////// 
double mth$dcos (double *) ; 

c_routine int math_cos (double *in, double *out) 
{ double a; 
real b; 

a = getdouble ( (real *) in) ; 
b = mth$dcos (&a) ; 
*{real *) out = b; 

return 1;} // math_cos 

/////////////////////////////////////////////////////////////////// 

// NAT INT TAN 

////////////////////////////////;////////////////////////////////// 
// Brief description: Returns a tanget of an angle. 
// Expected: in - passed real angle 
// Result: out - tanget of an angle 

///////////////////////// / f / s //////////////////////////////////// / / 
double mth$dtan (double *) ; 

coroutine void math_tan (real * in, real * out) 
{ double a, b; 
real r; 

a = getdouble (in) ? 
b = mth$dtan (&a) ; 
dtor (&b, &r) ; 

*out = r;} // math_tan 

/////////////////////////////////////////////////////////////////// 

>, NAT INT ASIN 

/////////////////////////////////////////////////////////////////// 
// Brief description: ASIN(x) returns the angle whose SIN is x. 
// Expected: in - passed real value 
// Result: out - angle 

//////////////////// / i f/ffffftffff/tt/ftf j i t i t i t i t i i t i i i i i i / i i i i i i / 
double mth$dasin (double *) ; 

coroutine void math_asin (double *in, double *out) 
{ double a; 
real b; 

a = getdouble ((real *) in) ; 
b = mth$dasin (&a) ; 

* (real *) out = b;} // math as in 

/////////////////////////////////////////////////////////////////// 

; NAT INT ACQS 

/////////////////////////// f ////////////////////////////////////// / 
fi Brief description: ACOS (x) returns the angle whose COS is x 
// Expected: in - passed real value 
// Result: out - angle . 

/////////////////////////////////////////////////////////////////// 
double mth$dacos (double *) / 

coroutine void math_acos (double *in, double *out) 
{ double a; 
real b; 

a = getdouble ((real *) in) ; 
b = mth$dacos (&a) ; 

*(real *) out = b;} // math acos 

/////////////////////////////////////////////////////////;///////// 
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« NAT INT ATAN 

/////////////////////////////////////////////////////////////////// 
// Brief description: ATAN (x) returns the angle whose TANGENT is x. 
// Expected: in - passed real value. 
// Result: out - angle. 

/////////////////////////////////////////////////////////////////// 
double mth$datan (double *) ; 

c_routine void math_atan (double *in, double *out) 
{ double a; 
real b; 

a = getdouble ((real *) in); 
b = mth$datan (&a) ; 

* (real *) out ■ b;} // math_atan 

/////////////////////////////////////////////// / ////////////////// / 

n NAT INT SINH 

t t ! / t i t t t J I t I I t I I 1 i I I t I I I I I I I I I I t I ! t I I I I I I I I I I 1 I I t I I i I I 1 I I t I 1 t 1 I I I t 

II Brief description: SINH(X) returns hyperbolic sine of X. 

// Expected: in - passed real value. 

// Result: out - hyperbolic sine of passed value 

/////////////////////////////////////////////////////////////////// 
double mth$dsinh (double *) ; 

coroutine void math_sinh (double *in, double *out) 
{ double a; 
real b; 

a = getdouble ((real *) in); 
b = mth$dsinh (&a) ; 

* {real *) out = b;} // math_sinh 

/////////////////////////////////////////////////////////////////// 

» NAT INT COSH ^ _ ~ 

//////////////////// 7 //////// t //////////// f //// t / t // t // t ///////// / i 
II Brief dexcription: 

// COSH(X) returns hyperbolic cosine of X. 
// Expected: in - passed real value. 

// Result: out - hyperbolic cosine of a passed value. 

/ / / / / / / / / / / / // // / / // / / // / / / / // // / / / / titi f i i i i ( / / ft ( t t /////// / j / t i / i 
double mth$dcosh (double *); 

coroutine void math_cosh (double *in, double *out) 
{ double a; 
real b; 

a = getdouble ( (real *) in) ; 
b = mth$dcosh (&a) ; 

* (real *) out « b?} // math_cosh 

I ) t t ) t t I t I f I I I I l i t I ( I t I I l I l I I I I I I I I I I 1 I I I I t I I I I ! I t t t I i I I I I t f I I ! f I t i 

n NAT INT TANH 

i t t t i i t t i i i t r i i t i t t t f t t t t t / t t / t t j t/// t / j i i t t i i i i i i i t i i t i t t i t i i i i i i i 

II Brief description: TANH (X) returns hyperbolic tangent of X. 
// Expected: in - passed real value. 

// Result: out - hyperbolic tanget of passed value. 

///////////////////////////////////^/////////////////////////////// 
double mth$dtanh (double *); 

coroutine void math_tanh (double *in, double *out) 
{ double a; 
real b; 

a = getdouble ( (real *) in) ; 
b = rath$dtanh (&a) ; 

Mreal *) out - b;} // math_tanh 

I i t I I I i f t i f f t t t t / t t t / t t t t t J J I 1 I I I t ! I t I i I t I i t I i t ! 1 I ! i t t ) I I I i f I I i I I I 1 

* NAT RND 

/////////////////////////////////////////////////////////////////// 
// Brief description: // Expected: value - numeric expression seed = global seed value 

// Result: Returns a random number between one and passed // numeric expression. 

///////////////////////////////////////y/////////////////;///////// 

void realip (real *, real *) ,- 

real mth$random (int *) ; 

void math_rnd (real * value) 

{ real f = mth$random (&seed) ; 

* value = f;} // math rnd 

i t i i i t t i t t f r t ( t t t f t t t t / / t t / t / t t i t i i t j i t t t t i t i t i t t t i i t i i i i i i i i i t i i i t 

n NAT SGN 

////////////////////////////////////y//////////////////y/////yy//// 

// Brief description: Returns the sign of a number 

// Expected: in - passed numeric expression 

// Result: out - sign of a a numeric expression 

// = +1 if expression > 0 

// = -1 if expression < 0 

// = 0 if expression = 0 

f t i / t f f f ( t f t f / / / / / / / / / t / i j / y / / y / / / y / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / y 
coroutine int mth$sgn (int *) ; 
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coroutine int realsgn (void *) ; 

void math_sgn (real * in, int *out) 
{ real n = 0; 

out[0j = $cmp$real (in, &n);} // math_sgn 

t i i ( i t t i t i i i t i t t t i i f i t t i i t i i i i i l i t t i i i i i t i i i t ( f r t f t f t t f / i f f / J J ) i t i J 



NAT INT ABS 



// y / y // y //// y /////// y //////// y //////////////////////////////////// ' 
// Brief expression: 

// Returns the absolute value of a numeric expression. 

// Expected: in - passed real expression 

// Result: out - the absolute value of real expression 

/////////////////////////////y///////////////////////////////////// 
double mth$dabs (double *) ; 
void math_abs (real * in, real * out) 
{ double a; 
real b,- 

a = getdouble ((real *) in); 
b = mth$dabs (&a) ,■ 

*(real *) out = b;} // math_abs 

/////////////////////////////////////////////////////////////////// 



m NAT REAL MAX 

i { t / / i i j j i i i i t i i t i i t i i t i i i i i i i i i i i t i i i i i i i t i t i i i i i i i i i i t i i i i i i i i i i i 
II Brief description: Returnes the larger of two reals X and Y. 
// Expected: x - real value 
// y - real value 

// Result: out - larger of two reals X and Y 

/////////////////////////////////////////////////////////////////// 
void math_real_max (real * x, real * y, real * out) 
{ if ($cmp$real (x, y) < 0) 
*out = *y; 
else 

*out - *x;} // math real max 

////////////// r ///////////////// / ///////////////////////////////// / 



NAT REAL MIN 

//////////////////////////////////////////////////// y / / / / 7 //////// / 



// Brief description: Returns the lesser of the two reals X and Y 
2? // Expected: x - passed real value // y passed real value // Result: out - lesser of two 

O reals X and Y 

s /////////////////////////////////////////////////////////////////// 

void math_real__min (real * x, real * y, real * out) 
s { if ($cmp$real (x, y) < 0) 
^ *OUt = *X; 

W else 

Ly *out = *y; } // math_real_min 

/////////// y / j ///////////// j ////////////////////////////////////// / 



t, NAT SQR 

t i //////////// t / /// f / / // f /// / ///////////// j //////// f ///////////// i / 
II Brief description: Returns a square root of a real value // Expected: in - passed real number 
// Result: out - square root 

/////////////////////////////////////////////////////////////////// 
double mth$dsqrt (double *) ; 
void mathjsqr (real * in, real * out) 
{ double a; 
real b; 

a = getdouble (in) ; 
b = mth$dsqrt (&a) ; 

*out = b;} // math_sqr 

I i f t / I I I f I I * i I I f I t I I i t l < I I t I t l I l I l t I I I t i ! I l f f ( f t / i / t / t t / ) f } / ) J f I ! I I 



NAT INT LOG 



/////////////////////////////////////////////////////////////////// 

II Brief description: Returns the natural logarithm of a real number // Expected: in - passed real value 

// Result: out- natural logarithm // 

/ / // // // // // // // // // // // // / / / // // // // / / / // // // // // // // // // // // // // j 
double mth$dlog (double *) ,- 
void math_log (real * in, real * out) 
{ double a; 
real b; 

a = getdouble (in); 
b = mth$dlog (&a) ? 

*out = b; } // math_log 

i { I f } I f > t J I l I I t i I I I I I I l I t i I ! I ! I I i I I f I I I I I I ) I I I 1 ( I I I I I I I I l i I f t t 1 I I t t 



// NAT LOG2 

/ / / / / // / / // / / / // / / // // // // // // // / / // // // i i j / // / / ///////// i //////// / 

// Brief description: Returns the base 2 logarithm of given real value 

// Expected: in - passed real value // Result: out - base 2 logarithm 

//////////////////////////////////////////////////////////////////y 

double mth$dlog2 (double *); 

void math_log2 (real * in, real * out) 
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{ double a; 
real b; 

a a getdouble (in) ; 
b = mth$dlog2 (&a) ; 

*out = b;} // math_log2 

//////////////////////////////;//////////////////////////////////// 

« NAT INT LOG10 

/////////////////////////////////////////////// f //////////////// s / ' 

(I Brief description: Returns the common logarithm of the real value 

// Expected; in - passed real value // Result: out - common logarithm 

/////////////////////////////////////////////////////////////////// 
double mth$dloglO (double *) ; 
void mathJLoglO (real * in, real * out) 
{ double a; 
real b; 

a = getdouble (in) ; 
b = mth$dloglO (&a) ; 

*out = b;} // math_loglO 
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